## Abstract

We present the energy minimization of atomic clusters as a promising problem class for continuous black box optimization benchmarks. Finding the arrangement of atoms that minimizes a given potential energy is a specific instance of the more general class of *geometry optimization* or *packing problems*, which are generally NP-complete. Atomic clusters are a well-studied subject in physics and chemistry. From the large set of available cluster optimization problems, we propose two specific instances: Cohn-Kumar clusters and Lennard-Jones clusters. The potential energies of these clusters are governed by distance-dependent pairwise interaction potentials. The resulting collection of landscapes is composed of smooth and rugged single-funnel topologies, as well as *tunable* double-funnel topologies. In addition, all problems possess a feature that is not covered by the synthetic functions in current black box optimization test suites: *isospectral symmetry*. This property implies that any atomic arrangement is uniquely defined by the pairwise distance spectrum, rather than the absolute atomic positions. We hence suggest that the presented problem instances should be included in black box optimization benchmark suites.

## 1 Introduction

For many years it has been a common shortcoming in the black box optimization community that the publication of novel search heuristics has often not been accompanied by rigorous benchmarks on balanced sets of test problems. The IEEE CEC competitions on real-valued black box optimization, as well as the GECCO workshops on black box optimization benchmarking (BBOB), constitute a welcome and much-needed effort that puts experimental research in search heuristics on more solid ground. Naturally, these test suites are composed of artificially designed test functions with certain mathematical properties such as (non-)separability, noise, scalability, multimodality or, in the case of quadratic forms, different condition numbers. Careful design of test functions allows the drawing of conclusions about certain invariance properties of the tested search algorithms. We argue that minimizing energy landscapes of atomic clusters provides a useful extension to the current CEC or BBOB test sets for two main reasons: First, these problems have the property of isospectral symmetry, a characteristic that is not covered by the current benchmark sets and that potentially affects different search heuristics in different ways. Second, atomic cluster problems can be considered real-world problems, since they model physical phenomena and share a similar problem structure with other important real-world optimization tasks such as, for example, sensor placement problems (Wu and Verma, 2008).

The key objective of the present article is to introduce black box optimization benchmark problems from energy minimization of atomic cluster configurations. The selection of problem instances is guided by two principles: First, the test suite should have benchmark problems of varying degrees of difficulty. Similar to synthetic benchmark suites that are composed of problems with unimodal or (highly) multimodal topology, we present benchmarks that can be seen as cluster analogs of these established functions. Second, all problem instances should be easy to implement, scalable to higher dimensions, and quick to evaluate on standard computers, thus allowing extensive benchmark runs with reasonable computational cost. To this end, we also make the MATLAB/Octave implementations of the present benchmarks publicly available on the authors’ website.

This article does not attempt any performance comparison between different optimization strategies. The included optimization runs solely serve as illustrations of how to use the benchmarks, and they shall provide a first reference for future comparisons of algorithms. We do, however, include a general introduction to atomic clusters and the associated minimization problems for computer scientists. Although most information can be found in standard physical chemistry textbooks, we feel that the provided information and references serve as a valuable entry point for researchers in evolutionary computation.

This article is structured as follows: In Section 2 we summarize the physical and computational aspects of atomic clusters. We also review standard *gradient-based* optimization approaches from chemistry for cluster problems. In Section 3 we outline general properties of the considered cluster instances and present a useful order parameter for cluster characterization. Section 4 introduces a novel cluster minimization problem based on Cohn-Kumar potentials and presents detailed information about its landscape features. Section 5 considers specific instances of Lennard-Jones clusters that are amenable to black box optimization. We conclude by summarizing the key results and providing ideas for further promising cluster instances.

## 2 Physical and Computational Aspects of Atomic Clusters

A cluster is a spatial arrangement of particles, typically a few tens to hundreds in number. In chemistry and physics, the study of clusters of atoms provides a means to understand nucleation phenomena. Nucleation describes the transition from a loose collection of atoms to the formation of bulk materials with particle numbers on the order of the Avogadro constant . Depending on chemical composition and physical conditions, different types of solids, such as crystals, quasi-crystals, amorphous solids, or glasses, emerge. The concept of energy landscapes has provided a fruitful view of these complex processes. This description considers the potential energy surface (PES) arising from the interacting atoms as a high-dimensional landscape with valleys, ridges, saddle points, and peaks (Stillinger and Weber, 1984; Wales and Scheraga, 1999). This concept is almost identical to Sewall Wright's fitness landscape metaphor (Wright, 1932), which has had tremendous impact in evolutionary biology and computation. It is now widely acknowledged that the global topology of a PES has an important influence on the way that atomic clusters and bulk material form. When the PES has a funnel-like shape where local minima are arranged in order of decreasing energy around the global minimum (or ground state), rapid evolution of the physical system toward this ground state is likely. Multi-funnel landscapes are composed of two or more such arrangements of local minima that are separated by large energy barriers. When a system is in one of the funnels, it is thus hard for it to escape and explore states in other funnels. Such landscapes are typical for glassy systems. Understanding and elucidating energy landscapes both in computer models and real experiments is thus an important research goal. In theoretical and computational approaches, PES can be discriminated by three fundamental model assumptions about the underlying physical system: (i) The number of particles in the system; (ii) the number of different atom types, and (iii) the classical or quantum-mechanical formulation of the energy. We first discuss the implications of these different model assumptions and state the choices we made in order to derive tractable benchmark problems for black box optimization.

### 2.1 Atomic Ensembles and Landscape Domains

The number of atoms in the system is the first important choice. Bulk systems with a huge number of atoms are approximated by a finite number of *N* atoms located in a rectangular prism (or a general parallelepiped) with periodic boundary conditions. In chemistry and crystallography, this rectangular prism is usually called the unit cell. When different kinds of atoms are considered, the unit cell should reflect the stoichiometry of the material, that is, the absolute ratio of the different atom types in the system. The correct size of the unit cell and the number of atoms it contains are usually not known a priori, thus hampering the setup of generic problem classes. The landscape domain of the described systems is determined by the number of atoms, the spatial extent of the unit cell, and a finite set of indicator variables that encode the atom types.

When the number of atoms is limited (usually to a few tens or hundreds), the spatial location of the atoms is bounded and no periodic boundary conditions are applied. This leads to atomic cluster or nano cluster systems. Over the past three decades, there has been increasing interest in such systems because of their wide-ranging chemical applications in catalysis, electronics, and energy conversion. We refer to Catlow et al. (2010) for an excellent chemical perspective on the topic. Many cluster systems have been studied, and the putative ground states (i.e., configurations of atoms attaining the global energy minimum) for different energy formulations are available, for instance in the Cambridge Cluster Database (CCD; Wales et al., 2009). Monatomic, binary, and ternary cluster systems are well studied, most prominently oxides such as zinc oxides or silica. Two different definitions of the space in which atomic clusters live are commonly used: Either the atoms populate the box-constrained Euclidean space, or they are restricted to the surface of the unit sphere where *d*_{E}(., .) denotes the Euclidean distance. For our benchmark problems, we consider both spaces. Finding particle distributions on the unit sphere that minimize a certain energy function is a long-standing problem since Thomson posed the question of how to optimally arrange electric charges on a sphere (Thomson, 1904). Optimal distributions of particles on the sphere are also known as spherical codes in coding theory. A collection of putative optimal spherical codes can be found in Sloane et al. (2000). Optimality is always defined with respect to a potential energy formulation.

### 2.2 Potential Energy Formulations and Ground States

*N*atoms in 3D space where the position of the

*i*th particle is denoted by

**p**

_{i}=(

*x*,

_{i}*y*,

_{i}*z*). Each configuration is restricted to a 3D box, that is, . The general form of an interatomic potential energy

_{i}*E*

_{IP}for many-body systems is: where represents the neighborhood of

**p**

_{i}. An established approach for the design of interatomic potentials is to split

*v*

_{IP}into two components. One component accounts for many-body energy contributions through pairwise interactions. The other component models the local environment of each individual atom. In general, the pairwise interaction term depends on all particles in the system; whereas the local environment is often defined within a specified distance range (neighborhood). Important examples of interatomic potentials that include both terms are the Finnis-Sinclair potential (Finnis and Sinclair, 1984) and its extension, the Sutton-Chen potential (Sutton and Chen, 1990), as well as tight-binding potentials (Cleri and Rosato, 1993). The famous Stillinger-Weber potential includes two-body and three-body interactions (Weber and Stillinger, 1985). The simplest instances of interatomic potentials are those that only consider pairwise interactions using isotropic pair potentials. The first potential of this kind dates back to John Lennard-Jones who introduced an empirical potential that describes the interaction between neutral atoms (Lennard-Jones, 1924). This Lennard-Jones (LJ) potential will be considered in Section 5. Another important instance is the Morse pair potential (Morse, 1929). Although all presented interaction potentials can be parameterized for different atom types, as well as atomic mixtures, we only consider monatomic clusters. The Thomson problem mentioned above can also be considered a minimization problem over a potential energy of pairwise interactions. There, the

*N*particles are electrons confined to the sphere that are arranged so as to minimize the total Coulomb potential.

The quality of potential energy models is usually assessed by their ability to either reproduce experimentally known ground states for various materials and clusters or to predict novel geometries as possible ground states that could guide experimentalists in their quest for novel forms of matter. Ground states define the macroscopic properties of the material. In a series of articles, Rechtsman, Stillinger, and Torquato (Rechtsman et al., 2005, 2006a, 2006b) introduced a new perspective on the topic: Instead of attempting to mimic nature as accurately as possible by tuning potential functions, they considered the inverse problem of how the shape of an isotropic pair potential must look in order to have as a ground state a predefined structure. In their simulation-guided optimization framework, they were able to design interaction potentials that result in bulk material with honeycomb (Rechtsman et al., 2005), square (Rechtsman et al., 2006a), cubic (Rechtsman et al., 2006b), diamond, and wurtzite (Rechtsman et al., 2007) lattices. All of their designed interaction potentials are based on multimodal isotropic pair potentials. This inverse statistical mechanics approach (Torquato, 2009) inspired Cohn and Kumar in their article, “Algorithmic design of self-assembling structures” (Cohn and Kumar, 2009) to derive convex pair potentials with provable ground states for clusters. We introduce two of these fundamental potentials in Section 4 and propose the resulting PES as a benchmark problem with a smooth single-funnel topology.

### 2.3 Standard Identification of Ground States

Atomic clusters and their ground states constitute an important problem class in chemistry and physics. From a computational complexity point of view, it is noteworthy that Wille and Vennik proved “that determining the ground state of a cluster of identical atoms, interacting under two-body central forces, belongs to the class of NP-hard problems” (Wille and Vennik, 1985, p. 419) using a proof by restriction. They showed that the cluster minimization problem includes the well-known traveling salesman problem as a special case. Adip (2005) revisited and refined the proof. Because of the existence of analytic gradients in most potential energy formulations, competitive algorithms to identify the ground states of larger clusters are hybrid stochastic-deterministic first-order methods that give no guarantees about the quality of the found solutions. Despite the tremendous number of different approaches to cluster optimization found in the (mostly chemistry or physics) literature, the basic ingredients of successful heuristics are few: a reasonably good global move set and an efficient quasi-Newton local gradient-based minimizer, mostly of the Broyden-Fletcher-Goldfarb-Shanno (BFGS) type. We refer to Wales's book (Wales, 2005) for further details, and to the excellent review by Hartke (2004) for an overview of applications of hybrid evolutionary algorithms. One of the most successful algorithms is the basin-hopping (BH) algorithm (Wales and Doye, 1997), which is based on (1) uniformly random variation of atomic positions, and (2) relaxation of the perturbed structure to the nearest local minimum using a gradient-based optimizer (either a conjugate gradient method, Wales and Doye, 1997, or BFGS, Wales, 2005). In this context, a basin is the collection of configurations that lead to the same local minimum for a given gradient-based minimization routine. Acceptance of a new structure is based on the Boltzmann criterion with respect to the previously accepted configuration. This is analogous to simulated annealing. In BH, however, the temperature parameter remains constant. Numerical results of BH and related memetic techniques for different cluster systems are scattered over hundreds of publications. For some systems, detailed knowledge about the number of minima, first-order saddle points, and global landscape structure is available. The present work is a first attempt to utilize this information for black box benchmarking.

## 3 Cluster Problems for Black Box Benchmarking

*n*of a cluster problem instance depends on the number of particles in the cluster and the space that the individual atoms populate. We consider problems not exceeding . Furthermore, we restrict ourselves to clusters of identical atoms with energy formulations based on scalar isotropic pair potentials

*u*

_{PP}. This leads to an energy formulation of the kind: where

*r*=

_{ij}*d*

_{E}(

**p**

_{i},

**p**

_{j}) is the Euclidean distance between particles

**p**

_{i}and

**p**

_{j}. The sole dependence of the isotropic pair potential on pairwise distances has several consequences. First, the resulting energy functions are nonseparable because all particles interact with all others. Second, they are, in general, continuously differentiable polynomial surfaces that are not bounded from above because most pair potentials diverge when

**p**

_{i}approaches any

**p**

_{j}. Third, the computational cost of energy function evaluation is low, but scales quadratically with the number of particles (and hence dimensions). Fourth, there is no unique set of absolute atomic positions that minimizes the energy; rather, all configurations that have a certain distance spectrum are energetically equivalent. This property is called isospectral symmetry and will be discussed in detail in Section 3.1. It is the key feature of cluster problems that is not represented in current standard benchmark suites.

A convenient property of atomic cluster problems is the intuitive illustration of (sub-)optimal solutions as a list of atomic coordinates that can be visualized in 3D space. In addition, the physical concept of order parameters enables low-dimensional descriptions of cluster configurations, as will be outlined in Section 3.2. Cluster problems also exhibit a wide variety of landscape topologies, thus leading to a rich and diverse benchmark set in the spirit of the CEC 2005 or BBOB test suites. Moreover, atomic cluster problems can be considered as real-world problems as they model physical systems. Other real-world problems, such as the sensor placement problem (Wu and Verma, 2008), share a similar problem structure. Algorithms that perform well on atomic cluster problems may therefore also prove successful for these applications.

### 3.1 Symmetry as a Novel Problem Characteristic

Two sources of symmetry arise in monatomic clusters: Consider the global minimum for a given energy function. When consists of the positions of *N* identical particles, *N*! possible permutations of atomic positions exist that attain the same ground state. Furthermore, the energy function does not discriminate between configurations that have identical sets of pairwise distances, that is, an identical distance spectrum. This characteristic, known as isospectral symmetry, implies that any transformation applied to that preserves all pairwise distances between the particles results in another minimum-energy configuration. A unique description of the cluster ground state is thus a set of *N*(*N*−1)/2 pairwise distances, rather than a set of coordinates **x**. Depending on the symmetry group of the minimum-energy configuration, different types of spectrum-preserving transformations exist. The simplest ones are translation, rotation, inversion, and reflection. However, it is also known that certain distance spectra can be generated by geometrically distinct configurations. Note that this type of symmetry is conceptually different from the symmetries found in standard benchmark functions, such as Rastrigin's function, where permutation symmetries often exist for local minima, but no symmetry exists for the global minimum. In the BBOB and CEC 2005 test suites, symmetries are also partially removed by (non-)linear transformations (Suganthan et al., 2005; Finck et al., 2009). Isospectral symmetry (and symmetry breaking) is a fundamental concept in physics, but its impact on black box optimization has been comparatively weak. We are only aware of the research of Van Hoyweghen and coworkers who analyzed symmetry “due to the interaction structure” (Van Hoyweghen and Naudts, 2000, p. 1073) in the problem. They describe the impact of such symmetry on the performance of evolutionary algorithms in the context of aggregated combinatorial problems (Van Hoyweghen and Naudts, 2000; Van Hoyweghen et al., 2002). They “claim that the occurrence of symmetry in the representation is another problem difficulty characteristic” (Van Hoyweghen et al., 2002, p. 317), and they propose different ways of modifying black box algorithms in order to cope with this difficulty. For continuous black box problems, such considerations are, to the best of our knowledge, so far missing. We hence propose the present benchmarks since it is conceivable that the performance of certain optimization algorithms (such as real-coded genetic algorithms) is more strongly impeded by isospectral symmetry than that of others (e.g., evolution strategies).

### 3.2 Bond Order Parameters for Cluster Characterization

**x**. Consider two configurations

**x**and

**y**, where

**y**is generated by the rotation of

**x**about the origin. Calculating

*d*

_{E}(

**x**,

**y**) would result in a value greater than 0, despite the fact that the two configurations are identical in terms of their pairwise distance spectrum. Measures that are invariant to isospectral symmetries provide a more robust way of characterizing the system. Designing such invariants for specific systems, however, is not trivial. Steinhardt and coworkers introduced a set of invariants for atomic cluster and bulk configurations that are applicable in our case:

*bond (orientational) order parameters*(BOP; Steinhardt et al., 1983). These parameters are indispensable for the analysis of nucleation phenomena and packing structure in bulk materials and cluster configurations. In this context, a bond does not refer to a covalent chemical bond, but is rather defined as the vector joining a pair of neighboring atoms. Neighborhood is defined by a distance cutoff. Bond-order parameters reflect the symmetry of bond orientations, regardless of absolute bond lengths. This is achieved by combination and normalization of certain spherical harmonics, resulting in the measures

*Q*and . The second-order invariants

_{l}*Q*are defined as: where and .

_{l}*N*denotes the number of bonds that are shorter than the cutoff distance

_{b}*r*

_{0}. The are spherical harmonics with being the polar and the azimuthal angle of the interatomic vector

**r**

_{ij}of length

*r*between atoms

_{ij}**p**

_{i}and

**p**

_{j}with respect to an arbitrary coordinate frame. The parameters are defined as: They are normalized versions of the third-order invariants where the coefficients are the so-called Wigner 3j symbols (Weisstein, 2010). Steinhardt and coworkers showed that the parameters

*Q*

_{4},

*Q*

_{6}, , and are sufficient for a detailed cluster shape spectroscopy (Steinhardt et al., 1983) of liquids, crystals, and glasses, since they discriminate between the most important symmetry groups. Specific atomic packings have unique combinations of values for this set. For instance, the parameter

*Q*4 discriminates between icosahedral (ico) and face-centered cubic octahedral (fcc) packing systems with values

*Q*

^{ico}

_{4}=0 and

*Q*

^{fcc}

_{4}=0.1909, respectively. In the context of black box optimization, we suggest using bond-order parameters as structural fingerprints of putative optimal cluster solutions, as well as for convenient visualization of optimization trajectories.

Alternative invariants can be constructed directly from the distance spectrum of the configurations. A wealth of such techniques exists in computer graphics and image processing for symmetry and shape descriptions. One instance is the concept of shape distributions (Osada et al., 2002), which might be interesting to examine in the context of black box optimization.

## 4 Cohn-Kumar Clusters

The first cluster instances we propose are based on pair potentials recently introduced by Cohn and Kumar (2009). Inspired by the inverse statistical mechanics approach, they designed pair potentials that result in provable ground states for certain cluster instances. To the best of our knowledge, there is no publication that considers these potentials. We hence introduce the name Cohn-Kumar (CK) potentials for the interaction potentials and Cohn-Kumar (CK) clusters for the resulting ground state clusters. The four pair potentials (CK1 – CK4) form the following provable minimum energy configurations:

An eight-particle CK1 cluster forms a 3D cube with six identical square faces.

A 20-particle CK2 cluster forms a 3D regular dodecahedron with 12 identical pentagonal faces.

A 16-particle CK3 cluster forms a 4D hypercube with eight identical cubic faces.

A 600-particle CK4 cluster forms a regular 120-cell in 4D with 120 dodecahedral faces.

*r*is the Euclidean distance between two particles. The CK2 pair potential is defined as: with

*t*=1−

*r*

^{2}/2. Both functions are designed to be monotonically decreasing and convex. Their graphs are depicted in Figure 1. For a system of

*N*particles with positions

**p**

_{i}, the energy functions are: and Again, and

*r*=

_{ij}*d*

_{E}(

**p**

_{i},

**p**

_{j}).

Let us first consider the two instances of CK1 and CK2 clusters for which the ground states have been proven by Cohn and Kumar (2009). The proofs are based on techniques from coding theory and linear programming. Details are given in Cohn and Kumar (2009). The optimal eight-particle CK1 (CK1_{8}) cluster configuration **x**^{CK1}_{min} is the 3D cube with identical square faces. This configuration belongs to the class of platonic solids and is depicted on the left of Figure 2. Due to the restriction of particle positions to the unit sphere, the following three unique distances occur: the cube edge length of , the face diagonal of length , and the cube diagonal length of 2. The energy minimum is . The dual polygon is the octahedron.

The optimal 20-particle CK2 (CK2_{20}) cluster configuration **x**^{CK2}_{min} is a dodecahedron with 12 identical pentagonal faces. Five unique distances occur in this platonic solid: The length *l*_{e} of the pentagonal edges is related to the radius *r*_{Sph} of the sphere on which the particles are located via . This leads to . The second, third, and fourth nearest neighbor distances have values of , , and , respectively. The largest occurring distance is the sphere diameter of 2. The minimum energy is . The dual polygon is the icosahedron. This dodecahedral configuration is depicted on the right of Figure 2.

Although these configurations are geometrically simple, it is hard to make them ground states of any classical pair potential energy function. Cohn and Kumar state the following (which is valid for all ground state configurations listed above): “The problem is that their facets are too large, which makes them highly unstable. Under ordinary potential functions, such as inverse power laws, these configurations are never even local minima, let alone global minima. In the case of the cube, one can typically improve it by rotating two opposite facets so they are no longer aligned. That lowers the energy, and indeed the global minimum appears to be the antiprism arrived at via a rotation (and subsequent adjustment of the edge lengths) (Cohn and Kumar, 2009, p. 9572). The square antiprism (or anticube) is, for instance, the minimum of the eight-electron Thomson problem, as well as the putative ground state for the eight-particle CK2 cluster. Despite the convexity of the CK pair potentials, the resulting energy functions are nonconvex even for the clusters with proven optimal configurations. Cohn and Kumar do not comment on the exact number of basins or the basin depths. For the optimal CK1 cluster they empirically found that out of 1,000 local optimization runs only six did not converge to the global minimum (Cohn and Kumar, 2009). Our own numerical results did not identify any local minima on the landscapes of the eight-atom CK1 cluster and the 20-atom CK2 cluster.

These results suggest two ways of using CK clusters as black box optimization benchmarks. The first way is solely based on the proven optimal configurations. A benchmark set with variable dimensionality can be constructed by fixing a number of optimal atomic positions, leading to well-defined multimodal problems. Consider the optimal eight-particle CK1 cluster. Fixing seven particles to optimal relative positions, we arrive at a 2D cluster problem with four basins of attraction: three that correspond to the square faces of equal size and one large basin consisting of the remaining three faces with the global minimum at the vacant corner position. Likewise, one could construct a 2D 10-funnel landscape from the optimal 20-particle CK2 cluster with nine suboptimal basins (the stable pentagons) and a large basin with the global minimum at the vacant position surrounded by three pentagonal faces. Decreasing the number of fixed atomic positions results in higher-dimensional problems with varying landscape structure.

The second way of constructing a standard benchmark set for black box optimization based on CK1 and CK2 clusters consists of varying the number of atoms on the sphere. We prefer this approach for its simplicity and suggest the following specification. We consider a CK1 and CK2 cluster with up to *N*=25 atoms. A system of *N* atoms confined to the surface of the unit sphere results in *n*=2*N* DOF. A natural landscape representation is based on spherical coordinates. The position **p**_{i} of atom *i* is defined by the pair of polar and azimuthal angles at unit radius, that is, . In order to construct a box-constrained problem, we restrict the angles to the interval , leading to . In principle, one could remove 4 DOF by fixing the position of a single atom on the sphere and fixing one polar and one azimuthal angle of one pair of atoms, thus removing certain symmetry properties of the problem. We do not, however, follow this approach here since we want to construct a benchmark with full isospectral symmetry. As in the CEC 2005 and BBOB test suite, we suggest the use of a budget of MAX_FES =10^{4}*n* for a single optimization run, and 25 independent runs per problem. We first focus on the eight-atom CK1 cluster and the 20-atom CK2 cluster, for which the optimal configurations and the corresponding energy values are known exactly. We then derive putative global optima for all other instances from numerical optimization runs in Section 4.2. These are used to specify a termination criterion for the level of solution accuracy . We also provide BOP values in order to characterize the symmetry of the optimal configurations.

### 4.1 Reference Black Box Experiments on the CK1_{8} and CK2_{20} Cluster Problems

In order to illustrate the use of black box optimizers on the eight-atom CK1 and the 20-atom CK2 cluster problems, and to provide a baseline for future performance comparisons, we run two sets of numerical optimization experiments. The first uses the covariance matrix adaptation evolution strategy (CMA-ES) with an increased population at each restart (IPOP-CMA-ES; Auger and Hansen, 2005). The second optimizer is MATLAB's function, implementing the Nelder-Mead (NM) simplex method (Nelder and Mead, 1965). CMA-ES is a stochastic optimization heuristic that iteratively samples a population of candidate solutions from a multivariate normal distribution. In each generation, CMA-ES adapts the mean and covariance matrix of the search distribution such as to increase the likelihood of sampling previously successful search directions. Additional adaptation of a step size parameter allows learning the appropriate scale of the problem. We refer to Hansen (2008) for further details. The NM simplex method is a deterministic direct-search heuristic that iteratively modifies the shape of an *n*-dimensional simplex so as to improve the worst candidate solution among the *n* simplex corners.

For IPOP-CMA-ES, we use the standard parameter settings reported by Auger and Hansen (2005; i.e., initial step size ), except for the parameter *incPop* that controls the population increase upon restart. We set this factor to 1.25 (instead of the standard value 2) since both our own experience and parameter tuning results (Smit and Eiben, 2010) showed superior performance when using a reduced factor on the CEC 2005 benchmarks as compared to the standard settings. For the NM algorithm, we perform independent restarts if it converges before reaching the global minimum, until the FES budget is exhausted. All initial configurations are drawn uniformly at random from the search domain. We use a solution accuracy of for both test problems (see also Table 5, discussed in Section 4.2).

The eight-atom CK1 cluster problem poses little challenge to either of the algorithms. Both methods almost always converge to the global minimum within the specified accuracy (CMA-ES with 100% and NM with 96% success rate). The statistics about the minimum energies reached by NM and CMA-ES are summarized in the left part of Table 1. The statistics about the number of FES that NM and CMA-ES needed to solve CK1_{8} are summarized in Table 2. From Table 2 we see that CMA-ES is more efficient in finding the ground state. While the minimum number of FES is comparable (5,317 vs. 5,879), the average number of FES (6,320.20 vs. 10,791.00) is considerably lower for CMA-ES.

. | CK1_{8}
. | CK2_{20}
. | ||
---|---|---|---|---|

Energy . | NM . | CMA-ES . | NM . | CMA-ES . |

Min | 6.343388 | 6.343388 | 746.730558 | 746.666667 |

7th | 6.343388 | 6.343388 | 746.775331 | 746.666667 |

Median | 6.343388 | 6.343388 | 746.825053 | 746.666667 |

19th | 6.343388 | 6.343388 | 746.872171 | 746.666667 |

Max | 6.343692 | 6.343388 | 746.910767 | 746.666667 |

Mean | 6.343400 | 6.343388 | 746.823867 | 746.666667 |

SD | 0.000061 | 0.000000 | 0.052601 | 0.000000 |

. | CK1_{8}
. | CK2_{20}
. | ||
---|---|---|---|---|

Energy . | NM . | CMA-ES . | NM . | CMA-ES . |

Min | 6.343388 | 6.343388 | 746.730558 | 746.666667 |

7th | 6.343388 | 6.343388 | 746.775331 | 746.666667 |

Median | 6.343388 | 6.343388 | 746.825053 | 746.666667 |

19th | 6.343388 | 6.343388 | 746.872171 | 746.666667 |

Max | 6.343692 | 6.343388 | 746.910767 | 746.666667 |

Mean | 6.343400 | 6.343388 | 746.823867 | 746.666667 |

SD | 0.000061 | 0.000000 | 0.052601 | 0.000000 |

Problem . | Min . | 7th . | Median . | 19th . | Max . | Mean . | SD
. | p_{s}
. |
---|---|---|---|---|---|---|---|---|

CK1_{8} | ||||||||

CMA-ES | 5,317 | 5,845 | 6,085 | 6,745 | 7,933 | 6,320.20 | 715.43 | 1 |

NM | 5,879 | 7,296 | 10,791 | 12,976 | — | 10,298.96 | 2,969.71 | .96 |

CK2_{20} | ||||||||

CMA-ES | 36,631 | 44,086 | 55,756 | 66,331 | 111,736 | 59,333.20 | 19,133.11 | 1 |

NM | — | — | — | — | — | — | — | — |

Problem . | Min . | 7th . | Median . | 19th . | Max . | Mean . | SD
. | p_{s}
. |
---|---|---|---|---|---|---|---|---|

CK1_{8} | ||||||||

CMA-ES | 5,317 | 5,845 | 6,085 | 6,745 | 7,933 | 6,320.20 | 715.43 | 1 |

NM | 5,879 | 7,296 | 10,791 | 12,976 | — | 10,298.96 | 2,969.71 | .96 |

CK2_{20} | ||||||||

CMA-ES | 36,631 | 44,086 | 55,756 | 66,331 | 111,736 | 59,333.20 | 19,133.11 | 1 |

NM | — | — | — | — | — | — | — | — |

The 20-atom CK2 cluster problem reveals a different picture: While CMA-ES always converges to the global minimum within the specified accuracy, Nelder-Mead fails in all 25 runs. The minimum energies reached by NM range from 746.730 to 746.910 (right part of Table 1). CMA-ES always reaches the perfect dodecahedral configuration with energy 746.666 without restarts.

Comprehensive structural information about the cluster configurations is provided by Steinhardt's bond-order parameters. In order to compute them, each cluster is augmented by a dummy atom placed at the origin. We suggest using a cutoff distance of *r*_{0}=1.01 for the calculation of BOP values in CK clusters. In Figure 3 we summarize the BOP values for the suboptimal configurations found by NM and the optimal ones reached by CMA-ES. We also depict the best configurations found by NM and CMA-ES. The BOP values for the suboptimal NM configurations indicate great structural diversity. Although the individual energy values are comparable, the BOP values vary considerably, especially for *Q*_{6} and . Among all BOP traces, the values of the best NM structure (highlighted in gray in Figure 3) are the closest ones to the BOP values of a perfect dodecahedron. This is confirmed by visual inspection of this structure, revealing a dodecahedron with distorted pentagonal faces (see right panel of Figure 3), as opposed to most other suboptimal structures that do not show any clear packing.

The BOP values can also be used to conveniently visualize the optimization path. In Figure 4 we show the trajectories of a typical CMA-ES run on the CK1_{8} and CK2_{20} cluster problems. For the CK1_{8} clusters, the *Q*_{4}, *Q*_{6}, , and values of CMA-ES’ mean converge to the optimal values after about 200 generations. For the CK2_{20} clusters, stable optimal BOP values for *Q*_{4}, *Q*_{6}, and are reached after about 2,600 generations. The values (data not shown) do not converge to the optimum CK2_{20} cluster. Such variation has also been observed for other instances. We therefore suggest to use the triplet *Q*_{4}, *Q*_{6}, and as a structural fingerprint for general CK clusters. For efficient restart strategies, it is conceivable to define (local) convergence of an algorithm in structure space in terms of these bond-order parameters rather than in terms of the original variables (i.e., the spherical coordinates).

The failure of the Nelder-Mead algorithm on the CK2_{20} cluster problem suggests that CK clusters are a nontrivial test case for black box optimizers. For now, we can only speculate about the poor performance of NM on this problem since NM and CMA-ES have the same invariance properties. Hansen tested a restart NM algorithm on the BBOB 2009 test suite and concluded that NM with restarts “allows searching unstructured multimodal landscapes comparatively effective, while a global topography within a multimodal or rugged landscape is not well exploited” (Hansen, 2009b). The robust performance of CMA-ES suggests a global single funnel topology of CK2_{20} with local minima on a smaller scale.

### 4.2 Putative Ground States of CK1 and CK2 Clusters for

So far we have analyzed two CK cluster instances with proven ground states. These instances have dimensionality *n*=16 and *n*=40, respectively. In order to construct a benchmark that spans a wider range of dimensions, we now consider all CK1 and CK2 clusters for even , that is, clusters containing up to *N*=25 atoms. Given the promising performance of CMA-ES on the previous instances, we use it as a tool for identifying putative ground states and low-energy local minima. We define the following computational experiments. For each instance, we run 25 standard CMA-ES runs without restarts until any of the standard convergence criteria are met. The initial step size is . The FES budget is restricted to MAX_FES = 10^{4}*n*. We store all putative global and local optima, energy values, and the number of FES CMA-ES needed to converge. We also calculate the values of the BOPs *Q*_{4}, *Q*_{6}, and for all observed structures.

We first report the results on the CK1 clusters. For these clusters, assuming an exponential scaling of the energies *E*_{CK1} of the putative ground states with cluster size yields the best least-squares fit (see Figure 5). The optimal fit is . The number of FES CMA-ES needed until convergence increases with *N*. The average number of FES scales linearly with cluster size. For , the variance is very low. For the other instances, however, some CMA-ES runs need considerably more FES to converge than others. This indicates that some problem instances exhibit considerably more complex landscapes than others, and that this phenomenon is not completely determined by the problem dimension. Nevertheless, all runs converge far before exhausting the FES budget. We summarize the information about putative CK1 ground states and low-energy minima as identified by CMA-ES in Table 3. We report the energies along with the BOP values for all minima. For *N*=10, 14, 16, 22, 23, 25, multiple low-energy minima were identified. From the wealth of generated data we discuss three instances in more detail.

N . | Energy . |
. Q_{4} |
. Q_{6} |
. |
---|---|---|---|---|

2 | 0.1084 | 1 | 1 | −0.0931 |

3 | 0.4630 | 0.3750 | 0.7408 | −0.0463 |

4 | 1.0583 | 0.5092 | 0.6285 | 0.0132 |

4 | 1.1029 | 0.5619 | 0.4369 | 0.0076 |

5 | 1.9838 | 0.6250 | 0.4556 | 0.0466 |

6 | 3.1501 | 0.7638 | 0.3536 | 0.0132 |

7 | 4.6241 | 0.5118 | 0.2861 | 0.0598 |

8 | 6.3434 | 0.5092 | 0.6285 | 0.0132 |

9 | 8.3580 | 0.1387 | 0.3561 | −0.0342 |

10A | 10.6645 | 0.2574 | 0.3289 | 0.0407 |

10B | 10.6646 | 0.2740 | 0.3651 | 0.0254 |

11 | 13.3828 | 0.0198 | 0.1129 | 0.1293 |

12 | 16.1847 | 0 | 0.0415 | −0.1698 |

13 | 20.5410 | 0.0736 | 0.2470 | −0.0032 |

14A | 25.1733 | 0.0771 | 0.2328 | 0.0131 |

14B | 25.1880 | 0.0729 | 0.2212 | 0.0111 |

14C | 25.2088 | 0.0278 | 0.0285 | −0.0931 |

15 | 30.9350 | 0.0332 | 0.1080 | 0.1196 |

16A | 38.4985 | 0.0464 | 0.0090 | 0.0931 |

16B | 38.5386 | 0.0615 | 0.2061 | −0.0631 |

17 | 47.3064 | 0.1044 | 0.0509 | 0.0931 |

18 | 59.9795 | 0.0052 | 0.1623 | −0.0931 |

19 | 78.2895 | 0.0487 | 0.1009 | 0.1476 |

20 | 94.1138 | 0.0630 | 0.1052 | −0.0407 |

21 | 122.3120 | 0.0173 | 0.1709 | −0.0274 |

22A | 151.7772 | 0.0294 | 0.0043 | 0.0132 |

22B | 153.1696 | 0.0307 | 0.0481 | 0.0026 |

23A | 202.9820 | 0.0139 | 0.1418 | −0.0329 |

23B | 203.0328 | 0.0145 | 0.1410 | −0.0277 |

24 | 236.1115 | 0.0164 | 0.0363 | 0.0132 |

25A | 314.0809 | 0.0047 | 0.1300 | −0.0158 |

25B | 314.0909 | 0.0048 | 0.1299 | −0.0176 |

25C | 321.7856 | 0.0235 | 0.1466 | −0.0316 |

N . | Energy . |
. Q_{4} |
. Q_{6} |
. |
---|---|---|---|---|

2 | 0.1084 | 1 | 1 | −0.0931 |

3 | 0.4630 | 0.3750 | 0.7408 | −0.0463 |

4 | 1.0583 | 0.5092 | 0.6285 | 0.0132 |

4 | 1.1029 | 0.5619 | 0.4369 | 0.0076 |

5 | 1.9838 | 0.6250 | 0.4556 | 0.0466 |

6 | 3.1501 | 0.7638 | 0.3536 | 0.0132 |

7 | 4.6241 | 0.5118 | 0.2861 | 0.0598 |

8 | 6.3434 | 0.5092 | 0.6285 | 0.0132 |

9 | 8.3580 | 0.1387 | 0.3561 | −0.0342 |

10A | 10.6645 | 0.2574 | 0.3289 | 0.0407 |

10B | 10.6646 | 0.2740 | 0.3651 | 0.0254 |

11 | 13.3828 | 0.0198 | 0.1129 | 0.1293 |

12 | 16.1847 | 0 | 0.0415 | −0.1698 |

13 | 20.5410 | 0.0736 | 0.2470 | −0.0032 |

14A | 25.1733 | 0.0771 | 0.2328 | 0.0131 |

14B | 25.1880 | 0.0729 | 0.2212 | 0.0111 |

14C | 25.2088 | 0.0278 | 0.0285 | −0.0931 |

15 | 30.9350 | 0.0332 | 0.1080 | 0.1196 |

16A | 38.4985 | 0.0464 | 0.0090 | 0.0931 |

16B | 38.5386 | 0.0615 | 0.2061 | −0.0631 |

17 | 47.3064 | 0.1044 | 0.0509 | 0.0931 |

18 | 59.9795 | 0.0052 | 0.1623 | −0.0931 |

19 | 78.2895 | 0.0487 | 0.1009 | 0.1476 |

20 | 94.1138 | 0.0630 | 0.1052 | −0.0407 |

21 | 122.3120 | 0.0173 | 0.1709 | −0.0274 |

22A | 151.7772 | 0.0294 | 0.0043 | 0.0132 |

22B | 153.1696 | 0.0307 | 0.0481 | 0.0026 |

23A | 202.9820 | 0.0139 | 0.1418 | −0.0329 |

23B | 203.0328 | 0.0145 | 0.1410 | −0.0277 |

24 | 236.1115 | 0.0164 | 0.0363 | 0.0132 |

25A | 314.0809 | 0.0047 | 0.1300 | −0.0158 |

25B | 314.0909 | 0.0048 | 0.1299 | −0.0176 |

25C | 321.7856 | 0.0235 | 0.1466 | −0.0316 |

The first instance is the CK1_{12} cluster. Its putative ground state is a Mackay icosahedron with 20 triangular faces (see Figure 12, discussed in section 5.1) with the well-known BOP pattern *Q*_{4}=0, *Q*_{6}=0.0415, and . The 13-atom Lennard-Jones cluster that is discussed in the next section exhibits the same symmetry with a central atom at the origin. From the CMA-ES runtimes, we see that CK1_{12} can be found rapidly (in less than FES on average) and robustly (all runs converge to the putative ground state).

For CK1_{14} clusters, CMA-ES converges to three different minima. In 21 out of the 25 runs, CMA-ES identifies the putative ground state (labeled 14A) with energy 25.1733. Two runs converge to a low-energy local minimum (labeled 14B) with energy 25.1880 and two runs to a local minimum (labeled 14C) with energy 25.2088. The corresponding configurations are shown in Figure 6. We highlight this cluster instance because the putative ground state might seem counterintuitive at first. A human observer would possibly favor structures 14B and 14C over 14A due to their apparent symmetry. They are, however, higher in energy than the putative ground state 14A. Moreover, structure 14A attains a value of , indicative of maximal cubic symmetry (Steinhardt et al., 1983).

The energy landscape of CK1_{16} clusters exhibits two competing low-energy structures as depicted in Figure 7. Eighteen out of the 25 CMA-ES runs converge to structure 16A, the remaining 7 runs find structure 16B. CK1_{16A} consists of two opposite rotated square faces (like the anticube), and triangular faces otherwise. CK1_{16B} has three square faces grouped around a central triangle, and otherwise triangular ones, leading to a different set of BOP values (see Table 3). In order to check whether the energy landscape explored by CMA-ES exhibits a single-funnel topology with two low energy minima at the bottom, or rather a double-funnel landscape, we conduct the following experiment. We start CMA-ES runs from the suboptimal low energy structure as the initial configuration with increasing initial values. These experiments reveal how much the initial configuration needs to be perturbed until CMA-ES is able to detect the putative globally optimal solution. We choose and repeat the experiment 50 times per . We monitor whether CMA-ES returns to the suboptimal solution or enters the putative ground state. The frequency of transition serves as an estimator for the transition probability P(CK1 CK1_{16A}) under CMA-ES exploration, and hence for the relative basin size of the suboptimal structure. The -dependent transition probability and a typical trajectory of CMA-ES in space leading from the suboptimal basin to the putative optimal basin are shown in Figure 8. The experiment suggests that below , CMA-ES does not leave the basin of the suboptimal solution. For larger the probability increases until it reaches a similar level as with global CMA-ES settings. Using , the probability is about 1/2 to fall into either minimum. The example CMA-ES trajectory for reveals an interesting pattern in space. Starting from the 16B structure, CMA-ES first performs a random walk until it clusters around configurations with BOP values , which most probably include the transition state between the two minima. The trajectory then smoothly converges to the final BOP values of the 16A structure. Trajectories that return to the suboptimal structure 16B behave similarly, with a cluster at , before smoothly converging back to the BOP values of the 16B configuration (data not shown). In summary, our experiments suggest that the CK1_{16} cluster landscape under CMA-ES exhibits a single-funnel topology with two competing minima at the bottom of the funnel. The putative optimal configuration is located in a considerably larger basin, hence representing a moderately difficult problem for CMA-ES. Nevertheless, it will be interesting to test other search heuristics on this problem, especially with respect to their susceptibility to the competing suboptimal solution.

We now present the key results for CK2 clusters. The scaling of the minimum energy with cluster size is shown in Figure 9. For the cluster sizes considered, the energies *E*_{CK2} of the putative ground states scale quadratically with cluster size. The best least-squares fit is achieved by . Note that this result can be explained by the moderate increase of the pair potential in the range of the observed distances for these cluster sizes (see Figure 1). The addition of a single particle hence only results in a quadratic increase of the energy. The average number of FES required by CMA-ES to converge to the putative minimum scales linearly with cluster size up to *N* = 12. For , the variance is very low. Compared to CK1 clusters, the number of FES needed to converge is considerably higher for larger CK2 clusters (well above FES on average), but is still an order of magnitude below the allowed FES budget. We summarize the information about putative ground states and local minima of CK2 clusters as identified by CMA-ES in Table 4. The energies and BOP values of all detected minima are reported. For , multiple minima were found. While these could be discussed analogously to our findings for CK1 clusters, we restrict ourselves to the interesting case of four-atom CK2 clusters. In 24 out of the 25 runs, CMA-ES finds the putative optimal ground state 4A, which is a regular tetrahedron with all pairwise distances equal to . In one run, CMA-ES finds a high-energy local minimum structure where the atoms form a pyramid with a larger triangular base face and three smaller triangular side faces (4B). The atoms form three distances of length and three distances of length . Both structures are depicted in Figure 10. It is surprising that the strictly convex CK2 pair potential produces a nonconvex energy landscape even in the four-atom case. Although the 4B structure is much higher in energy (3.0958) than the 4A tetrahedron (0.7901), its basin is relatively stable. Similar transition path experiments as were done for the CK1_{16} cluster reveal that when starting CMA-ES from the pyramidal structure, there is a high probability to converge back to the suboptimal structure even for as high as . Nevertheless, given the empirical hitting probability of 1/25, the overall basin size is negligible for CMA-ES.

V . | Energy . |
. Q_{4} |
. Q_{6} |
. |
---|---|---|---|---|

2 | 0 | 1 | 1 | −0.0931 |

3 | 0.0939 | 0.3750 | 0.7408 | −0.0463 |

4A | 0.7901 | 0.5092 | 0.6285 | 0.0132 |

4B | 3.0958 | 0.5312 | 0.5040 | 0.0048 |

5 | 6.0977 | 0.6250 | 0.4556 | 0.0466 |

6 | 12.0076 | 0.7638 | 0.3536 | 0.0132 |

7 | 29.2253 | 0.5536 | 0.0625 | −0.0931 |

8 | 49.3528 | 0.3736 | 0.2502 | −0.0931 |

9 | 75.8984 | 0.1502 | 0.3391 | −0.0436 |

10 | 108.7305 | 0.1702 | 0.1579 | −0.0931 |

11 | 149.0286 | 0.1043 | 0.3329 | −0.0515 |

12 | 192.0350 | 0 | 0.0415 | −0.1698 |

13A | 243.5499 | 0.0780 | 0.2679 | 0.0066 |

13B | 243.5531 | 0.1199 | 0.2714 | 0.0115 |

14 | 298.9590 | 0.2813 | 0.5036 | 0.0132 |

15 | 360.5035 | 0.0875 | 0.2977 | 0.0076 |

16A | 426.8723 | 0.0225 | 0.2791 | −0.0219 |

16B | 426.8726 | 0.0352 | 0.2890 | −0.0258 |

17 | 499.0473 | 0.1055 | 0.1994 | 0.0019 |

18 | 576.1469 | 0.1608 | 0.3495 | 0.0566 |

19A | 658.8684 | 0.1027 | 0.2900 | 0.1166 |

19B | 658.8689 | 0.0929 | 0.2822 | 0.1205 |

20 | 746.6667 | 0 | 0.2718 | 0.1698 |

21A | 840.1743 | 0.0389 | 0.1921 | 0.1141 |

21B | 840.1976 | 0.0372 | 0.1824 | 0.0835 |

21C | 840.2036 | 0.0364 | 0.1575 | 0.0922 |

22A | 938.8178 | 0.0353 | 0.1494 | 0.1517 |

22B | 938.8197 | 0.0326 | 0.1782 | 0.1098 |

23A | 1,042.8819 | 0.0335 | 0.0571 | −0.0206 |

23B | 1,042.8846 | 0.0193 | 0.0592 | −0.0507 |

23C | 1,042.8900 | 0.0185 | 0.1077 | −0.0846 |

23D | 1,042.9105 | 0.0181 | 0.1098 | 0.0301 |

24A | 1,152.1594 | 0.0058 | 0.0153 | 0.0132 |

24B | 1,152.1789 | 0.0037 | 0.0050 | −0.0931 |

25A | 1,266.8947 | 0.0130 | 0.0852 | 0.0706 |

25B | 1,266.9774 | 0.0186 | 0.0201 | 0.0931 |

V . | Energy . |
. Q_{4} |
. Q_{6} |
. |
---|---|---|---|---|

2 | 0 | 1 | 1 | −0.0931 |

3 | 0.0939 | 0.3750 | 0.7408 | −0.0463 |

4A | 0.7901 | 0.5092 | 0.6285 | 0.0132 |

4B | 3.0958 | 0.5312 | 0.5040 | 0.0048 |

5 | 6.0977 | 0.6250 | 0.4556 | 0.0466 |

6 | 12.0076 | 0.7638 | 0.3536 | 0.0132 |

7 | 29.2253 | 0.5536 | 0.0625 | −0.0931 |

8 | 49.3528 | 0.3736 | 0.2502 | −0.0931 |

9 | 75.8984 | 0.1502 | 0.3391 | −0.0436 |

10 | 108.7305 | 0.1702 | 0.1579 | −0.0931 |

11 | 149.0286 | 0.1043 | 0.3329 | −0.0515 |

12 | 192.0350 | 0 | 0.0415 | −0.1698 |

13A | 243.5499 | 0.0780 | 0.2679 | 0.0066 |

13B | 243.5531 | 0.1199 | 0.2714 | 0.0115 |

14 | 298.9590 | 0.2813 | 0.5036 | 0.0132 |

15 | 360.5035 | 0.0875 | 0.2977 | 0.0076 |

16A | 426.8723 | 0.0225 | 0.2791 | −0.0219 |

16B | 426.8726 | 0.0352 | 0.2890 | −0.0258 |

17 | 499.0473 | 0.1055 | 0.1994 | 0.0019 |

18 | 576.1469 | 0.1608 | 0.3495 | 0.0566 |

19A | 658.8684 | 0.1027 | 0.2900 | 0.1166 |

19B | 658.8689 | 0.0929 | 0.2822 | 0.1205 |

20 | 746.6667 | 0 | 0.2718 | 0.1698 |

21A | 840.1743 | 0.0389 | 0.1921 | 0.1141 |

21B | 840.1976 | 0.0372 | 0.1824 | 0.0835 |

21C | 840.2036 | 0.0364 | 0.1575 | 0.0922 |

22A | 938.8178 | 0.0353 | 0.1494 | 0.1517 |

22B | 938.8197 | 0.0326 | 0.1782 | 0.1098 |

23A | 1,042.8819 | 0.0335 | 0.0571 | −0.0206 |

23B | 1,042.8846 | 0.0193 | 0.0592 | −0.0507 |

23C | 1,042.8900 | 0.0185 | 0.1077 | −0.0846 |

23D | 1,042.9105 | 0.0181 | 0.1098 | 0.0301 |

24A | 1,152.1594 | 0.0058 | 0.0153 | 0.0132 |

24B | 1,152.1789 | 0.0037 | 0.0050 | −0.0931 |

25A | 1,266.8947 | 0.0130 | 0.0852 | 0.0706 |

25B | 1,266.9774 | 0.0186 | 0.0201 | 0.0931 |

In summary, we presented a detailed analysis of Cohn-Kumar clusters arising from two different strictly convex pair potentials. We analyzed the configurations where proven global minima exist and extended to other instances for up to *N*=25 atoms. In order to show the richness of the energy landscapes, we analyzed several instances in further detail using CMA-ES as a search heuristic and Steinhardt's bond-order parameters as measures to characterize the found energy minima. Table 5 summarizes our suggested CK cluster test suite settings.

## 5 Lennard-Jones Clusters

*r*=

_{ij}*d*

_{E}(

**p**

_{i},

**p**

_{j}) and

**p**

_{i}=(

*x*,

_{i}*y*,

_{i}*z*), the 3D position of atom

_{i}*i*. Parameter is the potential well depth (in units of energy) and is the equilibrium interatomic distance (in units of length) at zero temperature. Figure 11 depicts the unimodal shape of the LJ pair potential and the role of the parameters. The potential energy

*E*

_{LJ}of a cluster of

*N*LJ atoms is given by: Again, and

*r*=

_{ij}*d*

_{E}(

**p**

_{i},

**p**

_{j}). The ambient space of the atoms is the 3D Euclidean space. Knowledge about minimum-energy (or ground-state) configurations of LJ clusters allows predicting properties of crystallization or solid-liquid transitions of noble atomic mixtures at low temperatures. The presence of a well in the LJ pair potential implies that a collection of atoms faces the problem of frustration. Although all atoms “like” to have their neighbors at equilibrium distance, there is no geometric configuration that could achieve this for

*N*>4. This introduces multimodality in the landscape with competing low energy configurations.

For a long time it was believed that ground states of LJ clusters could be efficiently constructed by aufbau algorithms (Hoare, 1979). Starting from a seed structure with specific symmetry and *N*−1 atoms, these algorithms construct a putative ground state of *N* atoms by placing an additional atom at the energetically most favorable location and relaxing the resulting structure by energy-gradient descent. In the past 25 years, however, it has been shown that such algorithms are not able to identify many of today's known putative ground states. Wille (1987) identified the putative ground state of the LJ_{24} cluster using simulated annealing with a specialized problem-specific move set. This is the only ground state that has first been found by a method without the use of gradient information. Northby (1987) identified putative ground states for by searching lattices with icosahedral symmetry and gradient minimization. However, not all LJ instances follow an icosahedral symmetry, a fact that was only discovered in the mid-1990s. The ground states of LJ clusters with *N*=38, 75−77, 98, 101−103 atoms follow different packing schemes. They have mostly been identified by extensive application of unbiased gradient-based optimizers such as the hybrid genetic algorithm by Deaven and coworkers (Deaven et al., 1996) or basin-hopping (Li and Scheraga, 1987; Wales and Doye, 1997). The late discovery of these ground states is explained by the deceptive landscape topology with the global minimum located in a narrow funnel. In Section 5.2 we consider the archetypal (Wales and Scheraga, 1999) double funnel energy landscape of the LJ_{38} cluster and present a tuning technique that is able to smoothly deform the topology to a single funnel problem.

While large LJ clusters are prohibitive for black box optimization benchmarking due to their staggering number of local minima, we nonetheless propose small instances with up to atoms as meaningful benchmark problems. Minimizing the potential energy of a cluster of *N* atoms in 3D space defines a continuous optimization problem in *n*=3*N*−6 dimensions, since three translational and three rotational DOF can be removed from the system. This is achieved by placing the first atom at the origin of the Cartesian coordinate system, the second along the *x*-axis, and the third in the *x*−*y* plane. Hence, *N*=19 defines a problem in *n*=51 dimensions. We characterize the putative ground states of LJ clusters with in section 5.1 and present numerical optimization runs for selected instances. We emphasize that LJ clusters have, despite their widespread use in chemistry and physics, not been subject to rigorous studies using black box heuristics. We are only aware of two publications where certain small cluster instances have been optimized using evolutionary algorithms: Müller and coworkers presented some initial results for CMA-ES on LJ clusters with *N*=8, 27 (Müller et al., 2003); and Call and coworkers optimized LJ_{26} with a specialized PSO (Call et al., 2007).

### 5.1 Lennard-Jones Clusters for

All known putative ground states of LJ clusters for the standard parameterization are available from the Cambridge Cluster Database (Wales et al., 2009). We characterize the structures using the BOP parameters *Q*_{4}, *Q*_{6}, and with the suggested *r*_{0}=1.391 for LJ clusters (Doye et al., 1999b). The data are summarized in Table 6. Most of the structures have a low *Q*_{4} value, indicating icosahedral symmetry (Steinhardt et al., 1983). As previously mentioned, the LJ_{13} ground state is identical (with identical BOP values) to the putative CK1_{12} ground state with an additional central atom at the origin (see also Figure 12). This instance is also the first of the so-called magic number structures. Magic number LJ clusters are based on complete multilayer Mackay icosahedra with atoms. They are so stable that they are regularly found in NMR experiments, for instance, of Xeon clusters (Hasse, 1991). The number of energy minima in these clusters scales exponentially with cluster size. For the largest cluster considered here, the LJ_{19} instance, Kunz and Berry (1995) estimated the number of geometrically distinct local minima to be around half a million.

N
. | Energy . |
. Q_{4} |
. Q_{6} |
. |
---|---|---|---|---|

3 | −3 | 0.3750 | 0.7408 | −0.0463 |

4 | −6 | 0.1909 | 0.5745 | −0.0132 |

5 | −9.1039 | 0.0013 | 0.4297 | 0.0314 |

6 | −12.7121 | 0.1909 | 0.5745 | −0.0132 |

7 | −16.5054 | 0.0148 | 0.1604 | −0.0931 |

8 | −19.8215 | 0.0644 | 0.1467 | 0.0015 |

9 | −24.1134 | 0.0156 | 0.1226 | −0.0391 |

10 | −28.4225 | 0.0382 | 0.1359 | −0.0349 |

11 | −32.7660 | 0.0298 | 0.1384 | 0.0485 |

12 | −37.9676 | 0.0178 | 0.1186 | 0.1209 |

13 | −44.3268 | 0 | 0.0415 | −0.1698 |

14 | −47.8452 | 0.0283 | 0.0437 | 0.0128 |

15 | −52.3226 | 0.0069 | 0.0437 | −0.0790 |

16 | −56.8157 | 0.0208 | 0.0653 | −0.0888 |

17 | −61.3180 | 0.0216 | 0.0849 | −0.0778 |

18 | −66.5309 | 0.0148 | 0.0676 | 0.1613 |

19 | −72.6598 | 0.0043 | 0.0056 | 0.0931 |

N
. | Energy . |
. Q_{4} |
. Q_{6} |
. |
---|---|---|---|---|

3 | −3 | 0.3750 | 0.7408 | −0.0463 |

4 | −6 | 0.1909 | 0.5745 | −0.0132 |

5 | −9.1039 | 0.0013 | 0.4297 | 0.0314 |

6 | −12.7121 | 0.1909 | 0.5745 | −0.0132 |

7 | −16.5054 | 0.0148 | 0.1604 | −0.0931 |

8 | −19.8215 | 0.0644 | 0.1467 | 0.0015 |

9 | −24.1134 | 0.0156 | 0.1226 | −0.0391 |

10 | −28.4225 | 0.0382 | 0.1359 | −0.0349 |

11 | −32.7660 | 0.0298 | 0.1384 | 0.0485 |

12 | −37.9676 | 0.0178 | 0.1186 | 0.1209 |

13 | −44.3268 | 0 | 0.0415 | −0.1698 |

14 | −47.8452 | 0.0283 | 0.0437 | 0.0128 |

15 | −52.3226 | 0.0069 | 0.0437 | −0.0790 |

16 | −56.8157 | 0.0208 | 0.0653 | −0.0888 |

17 | −61.3180 | 0.0216 | 0.0849 | −0.0778 |

18 | −66.5309 | 0.0148 | 0.0676 | 0.1613 |

19 | −72.6598 | 0.0043 | 0.0056 | 0.0931 |

We now provide numerical experiments on three select instances: LJ_{7}, LJ_{13}, and LJ_{19}. The putative ground state structures of these clusters are depicted in Figure 12. We focus on these three instances since they have been extensively analyzed in Wales (2005). LJ_{13} and LJ_{19} are also considered in Doye et al. (1999a). Detailed information is available about the number of unique local minima, first-order and higher-order saddle points, as well as disconnectivity graphs of low-energy minima. Prior analysis also showed that the energy landscapes of these clusters exhibit single funnel topologies, which should make them feasible for many black box optimizers, despite the huge number of local minima (Doye et al., 1999a).

We imagine these instances as cluster analogs of Rastrigin-like test functions. For the numerical experiments, we suggest the benchmark settings summarized in Table 8. We illustrate the performance of a black box optimizer on these test cases by running IPOP-CMA-ES with standard settings ( in this benchmark case), yet again with a smaller increase factor for the population size (*incPop*=1.25; restarts occurred in all LJ cluster test cases). The statistics about the number of FES and the success rate *p*_{s} for CMA-ES runs on LJ_{7}, LJ_{13}, and LJ_{19} are summarized in Table 7. CMA-ES always finds the optimal structure of LJ_{7}. In 24 out of the 25 runs, it also solves the LJ_{13} cluster problem. Inspection of Doye's and Wales's disconnectivity graphs for LJ_{13} (Doye et al., 1999a) reveals that CMA-ES converges to the second-best optimum in the one case where it does not find the global optimum. The LJ_{19} problem is more challenging. Only in two out of the 25 runs is CMA-ES able to identify the putative ground state. In the other runs, CMA-ES detects twice the second-best minimum, and in the remaining cases, a diverse set of local low energy minima. These results confirm that LJ clusters with moderate numbers of atoms are amenable for black box optimization, but are considerably more difficult for CMA-ES than CK1/CK2 clusters with the same number of degrees of freedom. In benchmark scenarios with a limited function evaluation budget, LJ clusters with larger numbers of atoms are likely to be hard for black box search.

. | Min . | 7th . | Median . | 19th . | Max . | Mean . | std . | p_{s}
. |
---|---|---|---|---|---|---|---|---|

LJ_{7} | 3,805 | 5,509 | 15,108 | 23,980 | 89,819 | 18,220.36 | 18,218.27 | 1 |

LJ_{13} | 13,007 | 51,827 | 101,726 | 178,359 | — | 109,377.12 | 80,713.15 | 0.96 |

LJ_{19} | 31,726 | — | — | — | — | 270,870.50 | 338,201.40 | 0.08 |

. | Min . | 7th . | Median . | 19th . | Max . | Mean . | std . | p_{s}
. |
---|---|---|---|---|---|---|---|---|

LJ_{7} | 3,805 | 5,509 | 15,108 | 23,980 | 89,819 | 18,220.36 | 18,218.27 | 1 |

LJ_{13} | 13,007 | 51,827 | 101,726 | 178,359 | — | 109,377.12 | 80,713.15 | 0.96 |

LJ_{19} | 31,726 | — | — | — | — | 270,870.50 | 338,201.40 | 0.08 |

### 5.2 The LJ_{38} Cluster as a High-Dimensional Benchmark with Tunable Landscape Topology

_{38}cluster problem, which has also been the inspiration for Lunacek's double funnel test case (Whitley, 2010). We provide a complete description of the proposed benchmark. An initial description with extensive numerical optimization runs using IPOP-CMA-ES was reported by Müller and Sbalzarini (2009). The topology of the standard LJ

_{38}energy landscape was widely studied in the literature (Barron et al., 1996; Doye et al., 1999b; Wales, 2004). It exhibits a double funnel structure where the global minimum-energy configuration with fcc octahedral symmetry lies in the narrow funnel. The majority of local minima, most of them with icosahedral (ico) symmetry, populate the wider funnel. Figure 13a shows a sketch of this landscape. Leary (2000) estimated the size of the optimal funnel from monotonic sequence BH runs to be around 12.4% of the entire configuration space. For 1,000 randomly generated configurations Leary applied the BH move set and only accepted improving configurations. In 124 out of the 1,000 runs, BH reached the fcc structure. For standard BH, however, the probability is much lower, since most runs converge to structures within the larger, suboptimal funnel. In fact, the optimal fcc structure was not found by unbiased optimization methods, but by a combination of geometric intuition and local minimization (Barron et al., 1996). Barron and coworkers first constructed initial LJ

_{38}configurations on fcc lattices and then applied gradient-based minimization in order to relax these structures under the LJ pair potential (Barron et al., 1996). Doye and coworkers (Doye et al., 1999b) eventually revealed the paradigmatic double-funnel nature of the landscape. They characterized it in terms of the number and location of minima, structural diversity of the minima, and the energy barrier between the two funnels. Based on this information, Doye realized that the problem of finding the global minimum of LJ

_{38}can be simplified by introducing a penalty term that simulates compression in the original energy function (Doye, 2000). Compression can be seen as a transformation of the PES that favors more compact structures. This broadens the funnel that contains the more compact fcc structures and narrows the ico funnel. Doye proposed the following penalized energy function: where and where

*p*

_{cm}is the center of mass of the cluster. The compression term has the form of an atomic positional variance. For the best icosahedral structure, it is

*Q*

_{c}(

**x**

^{ico})=96.1624; and for the best fcc structure

*Q*

_{c}(

**x**

^{fcc})=91.6369, hence distorting the energy difference between the two competing structures in favor of the octahedron. The scalar parameter controls the magnitude of compression. The effect of the -dependent compression on the topology of the PES is visualized with disconnectivity graphs in Doye (2000) and Wales (2004, pp. 338–339). When , we recover the original LJ

_{38}cluster problem. For , the PES exhibits a clear single-funnel topology. We sketch this phenomenon in Figure 13. Although the compression term also lowers the average barrier between local minima (Doye, 2000), the system still contains a staggering number of local minima for all considered values.

We propose the 38-atom LJ cluster with the energy function defined in Equation (13) as a high-dimensional tunable test case to study the performance of black box optimization methods as a function of landscape topology. In particular, we suggest the following benchmark scenario. First, the budget of allowed FES should be considerably increased. Second, an algorithm should be tested with . The lower the for which the algorithm can still find the putative ground state, the less sensitive it is to the landscape topology. We obtain the putative ground states for the different using reference CMA-ES runs. Therefore, we start CMA-ES from the known optimal **x**^{fcc} structure of the original problem with a small . From there, it quickly converges to slightly different fcc structures for varying . The energies and BOP values of the putative ground states are summarized in Table 9. For , the resulting energies match the ones reported in Doye (2000), providing further confidence that the ground states found here are correct. Visual inspection and the computed BOP values indicate that the structures minimizing the modified energy function in Equation (13) are almost identical to the optimal fcc configuration of the original, uncompressed problem. Doye and coworkers (Doye et al., 1999b) show that, among the different BOPs, *Q*_{4} is best suited for discriminating between fcc and ico structures, with *Q*_{4}(**x**^{ico})=0 and *Q*_{4}(**x**^{fcc})=0.1909. Table 10 summarizes our proposed specification for the tunable LJ_{38} cluster benchmark. The bounds are far from being tight with respect to the optimal structure. Both **x**^{ico} and **x**^{fcc} would fit into the [−2, 2]^{n} box. We propose the larger bounds for two reasons: First, we want to minimize effects from boundary handling techniques. Second, we want to test the capability of the optimization algorithm to cope with uninformative regions of parameter space. Enlarging the box adds plateau-like regions to the energy landscape, because a particle that is far away from the cluster experiences only a small force that draws it toward the cluster.

. | Energy . |
. Q_{4} |
. Q_{6} |
. |
---|---|---|---|---|

0 | −173.92843 | 0.19090 | 0.57446 | −0.01316 |

0.5 | −128.42914 | 0.19090 | 0.57445 | −0.01316 |

1 | −83.50595 | 0.19090 | 0.57444 | −0.01316 |

1.5 | −39.08437 | 0.19091 | 0.57443 | −0.01316 |

2 | 4.89246 | 0.19091 | 0.57442 | −0.01316 |

2.5 | 48.46946 | 0.19091 | 0.57441 | −0.01316 |

3 | 91.68310 | 0.19092 | 0.57440 | −0.01316 |

3.5 | 134.56362 | 0.19092 | 0.57439 | −0.01316 |

4 | 177.13651 | 0.19092 | 0.57438 | −0.01316 |

4.5 | 219.42358 | 0.19092 | 0.57437 | −0.01316 |

5 | 261.44371 | 0.19092 | 0.57436 | −0.01316 |

. | Energy . |
. Q_{4} |
. Q_{6} |
. |
---|---|---|---|---|

0 | −173.92843 | 0.19090 | 0.57446 | −0.01316 |

0.5 | −128.42914 | 0.19090 | 0.57445 | −0.01316 |

1 | −83.50595 | 0.19090 | 0.57444 | −0.01316 |

1.5 | −39.08437 | 0.19091 | 0.57443 | −0.01316 |

2 | 4.89246 | 0.19091 | 0.57442 | −0.01316 |

2.5 | 48.46946 | 0.19091 | 0.57441 | −0.01316 |

3 | 91.68310 | 0.19092 | 0.57440 | −0.01316 |

3.5 | 134.56362 | 0.19092 | 0.57439 | −0.01316 |

4 | 177.13651 | 0.19092 | 0.57438 | −0.01316 |

4.5 | 219.42358 | 0.19092 | 0.57437 | −0.01316 |

5 | 261.44371 | 0.19092 | 0.57436 | −0.01316 |

Problem . | LJ_{38} with compression
. |
---|---|

Runs/problem | 25 |

n | 108 |

MAX_FES | |

Termination | If FES=MAX_FES or |

Initialization and bounds | Uniformly at random in [−4, 4]^{n} |

Problem . | LJ_{38} with compression
. |
---|---|

Runs/problem | 25 |

n | 108 |

MAX_FES | |

Termination | If FES=MAX_FES or |

Initialization and bounds | Uniformly at random in [−4, 4]^{n} |

Numerical experiments on this benchmark show that IPOP-CMA-ES can solve it for , but not for , even when using billions of function evaluations (Müller and Sbalzarini, 2009). It is an open question whether CMA-ES variants that are designed for multi-funnel problems, such as the particle swarm CMA-ES (Müller et al., 2009) or BIPOP-CMA-ES (Hansen, 2009a), can cope with this problem. To date, no gradient-free black box optimization algorithm has been reported to solve the LJ_{38} test case without compression. Since many real-world applications involve multi-funnel landscapes, we believe that the tunable LJ_{38} problem with varying degree of compression presents a challenging test case for the black box optimization community that might prove instrumental in the design and analysis of new search heuristics.

## 6 Discussion and Conclusions

Finding the optimal spatial arrangement of atoms that minimizes the potential energy of a cluster system constitutes a promising problem class for continuous black box optimization benchmarks. To this end, we presented several atomic cluster problems and analyzed their energy landscapes and putative optima. We specifically proposed Cohn-Kumar (CK) clusters and Lennard-Jones (LJ) clusters, whose energies are given by sums over distance-dependent pairwise potentials.

We showed that CK clusters exhibit smooth landscapes with a single global minimum or few local minima. On CK cluster instances for which proven ground states are known, we illustrated the performance of a restart Nelder-Mead simplex algorithm and of IPOP-CMA-ES. IPOP-CMA-ES outperformed Nelder-Mead in terms of robustness, speed, and solution quality. For all other CK clusters up to *N*=25, we found putative global minima and several low-energy local minima from extensive numerical simulations. This provides the necessary information for a benchmark suite in up to *n*=50 dimensions.

The presented LJ cluster instances are known to exhibit rugged single-funnel topologies, as well as tunable double-funnel topologies. IPOP-CMA-ES was able to identify putative ground states for LJ clusters up to *N*=19. We further proposed the 38-atom LJ cluster with compression as a benchmark to assess the sensitivity of optimization heuristics with respect to landscape topology.

All cluster problems possess isospectral symmetry as a novel characteristic that is not covered by the test functions in current black box benchmark suites. They hence allow for an assessment of whether and how well black box algorithms can cope with this problem feature. We suggested using bond-order parameters as symmetry-invariant measures to characterize and compare structures. Search trajectories of black box optimizers can also be conveniently represented in the space spanned by these parameters.

We note that the present benchmarks could be extended to be composed of further cluster problems with properties that are not represented in the proposed benchmark set. We believe that two meaningful extensions would be Morse clusters and minimum second-moment sphere packings. Morse clusters can be used to design fixed dimensional problems with a tunable degree of ruggedness, but identical global topology. For the 13-atom Morse cluster, a thorough characterization of the landscape topography in function of a single parameter of the Morse pair potential can be found in Cox et al. (2006). Putative ground states of Morse clusters for and different values are reported in the Cambridge Cluster Database (Wales et al., 2009). The minimum second-moment sphere packing problem consists of finding the configuration of nonoverlapping spheres of identical radius that has the smallest second moment of the positions about the center of mass (Sloane et al., 1995). This defines a convex objective function with nonconvex quadratic constraints. The objective function is, in fact, identical to the compression penalty *Q*_{c} of the LJ_{38} test case. The hard-sphere constraints, however, turn the minimum second-moment sphere packing problem into a discontinuous problem where derivatives do not exist. Sloane and coworkers applied a variety of methods to this problem, ranging from direct search heuristics (simulated annealing) to complete enumeration. In the original article, they presented putative optimal configurations up to *N*=32. Putative optimal configurations up to *N*=99 are listed, elsewhere (Sloane et al., 1997). Both theoretical and experimental results (Meng et al., 2010) suggest that the energy landscapes of sphere-packing problems are not strongly funneled, but contain distinct local minima that are separated by large barriers. We expect that confirming or improving the currently known putatively optimal finite sphere packings is a formidable challenge for black box optimization methods.

In summary, we propose that atomic cluster problems should be included in black-box benchmark suites in order to better assess the efficacy, efficiency, and generality of modern search heuristics. We invite all researchers in the black-box optimization community to test their favorite algorithms on the presented test cases. Along with this article we will release publicly available MATLAB/Octave code that implements the presented cluster problems and the relevant bond-order parameter calculations.^{1}

## Acknowledgments

We thank the three anonymous reviewers for their comments which considerably improved the quality of the manuscript.

## References

*N*-body potential for transition metals

*GECCO ’09*

*j*-symbol. From MathWorld—A Wolfram Web Resource