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 dE(., .) 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

When a collection of particles interacts on the atomic scale, the potential energy landscape arises from the forces between all electrons and nuclei of the atoms. In order to arrive at a tractable model of the PES, the fundamental assumption in both classical and quantum mechanical (QM) formulations is the Born-Oppenheimer approximation (Born and Oppenheimer, 1927): Based on the large discrepancy between nucleus and electron masses, this approximation allows for a separation of the energy into an electronic and a nuclear component. Existing formulations of potential energy differ in the way they describe these components. Depending on the specific properties of the atoms in the system, different levels of detail are necessary in realistic models. Two main models gained popularity in the past decades: QM electronic structure techniques and classical interatomic potentials. Despite today's increasing availability of high-performance computing environments, QM-based energy calculations of atomic ensembles are still very costly. Classical interatomic potentials offer a convenient alternative to calculate the potential energy of many-particle systems in a rapid manner. The principle underlying the design of these potentials is largely empirical. One attempts to reproduce the experimentally observed dynamics and ground states of specific types of matter through careful parameterization of simple analytical functions. Consider a cluster of N atoms in 3D space where the position of the ith particle is denoted by pi=(xi, yi, zi). Each configuration is restricted to a 3D box, that is, . The general form of an interatomic potential energy EIP for many-body systems is:
formula
1
where represents the neighborhood of pi. An established approach for the design of interatomic potentials is to split vIP 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

The physical and computational foundations of atomic clusters can be exploited for the development of a cluster benchmark library. We adopt several design criteria of the established CEC 2005 (Suganthan et al., 2005) and BBOB test suites (Finck et al., 2009). The dimensionality 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 uPP. This leads to an energy formulation of the kind:
formula
2
where rij=dE(pi, pj) is the Euclidean distance between particles pi and pj. 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 pi approaches any pj. 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

Due to the presence of isospectral symmetry, it can be misleading to describe and compare cluster configurations in absolute coordinates x. Consider two configurations x and y, where y is generated by the rotation of x about the origin. Calculating dE(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 Ql and . The second-order invariants Ql are defined as:
formula
3
where
formula
4
and . Nb denotes the number of bonds that are shorter than the cutoff distance r0. The are spherical harmonics with being the polar and the azimuthal angle of the interatomic vector rij of length rij between atoms pi and pj with respect to an arbitrary coordinate frame. The parameters are defined as:
formula
5
They are normalized versions of the third-order invariants
formula
6
where the coefficients are the so-called Wigner 3j symbols (Weisstein, 2010). Steinhardt and coworkers showed that the parameters Q4, Q6, , 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 Q4 discriminates between icosahedral (ico) and face-centered cubic octahedral (fcc) packing systems with values Qico4=0 and Qfcc4=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:

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

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

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

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

We restrict ourselves to the first two potentials with ground state clusters living in 3D space. In both cases, the ambient space of the particles is , leading to 2 DOF per particle. The CK1 pair potential is defined as:
formula
7
where r is the Euclidean distance between two particles. The CK2 pair potential is defined as:
formula
8
with t=1−r2/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 pi, the energy functions are:
formula
9
and
formula
10
Again, and rij=dE(pi, pj).
Figure 1:

The Cohn-Kumar (CK) potentials uCK1(r) and uCK2(r) in log scale versus distance r for . In linear scale both pair potentials are strictly convex.

Figure 1:

The Cohn-Kumar (CK) potentials uCK1(r) and uCK2(r) in log scale versus distance r for . In linear scale both pair potentials are strictly convex.

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 (CK18) cluster configuration xCK1min 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.

Figure 2:

Configurations of the optimal eight-atom CK1 cluster (left) and the optimal 20-atom CK2 cluster (right).

Figure 2:

Configurations of the optimal eight-atom CK1 cluster (left) and the optimal 20-atom CK2 cluster (right).

The optimal 20-particle CK2 (CK220) cluster configuration xCK2min is a dodecahedron with 12 identical pentagonal faces. Five unique distances occur in this platonic solid: The length le of the pentagonal edges is related to the radius rSph 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=2N DOF. A natural landscape representation is based on spherical coordinates. The position pi 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 =104n 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 CK18 and CK220 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 CK18 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.

Table 1:
Statistics about the minimum energies reached by NM and CMA-ES for the CK18 and CK220 cluster problems.
CK18CK220
EnergyNMCMA-ESNMCMA-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 
CK18CK220
EnergyNMCMA-ESNMCMA-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 
Table 2:
Statistics about the number of FES used by NM and CMA-ES to reach the optimal configurations for the CK18 and CK220 cluster problems.
ProblemMin7thMedian19thMaxMeanSDps
CK18         
   CMA-ES 5,317 5,845 6,085 6,745 7,933 6,320.20 715.43 
   NM 5,879 7,296 10,791 12,976 — 10,298.96 2,969.71 .96 
CK220         
   CMA-ES 36,631 44,086 55,756 66,331 111,736 59,333.20 19,133.11 
   NM — — — — — — — — 
ProblemMin7thMedian19thMaxMeanSDps
CK18         
   CMA-ES 5,317 5,845 6,085 6,745 7,933 6,320.20 715.43 
   NM 5,879 7,296 10,791 12,976 — 10,298.96 2,969.71 .96 
CK220         
   CMA-ES 36,631 44,086 55,756 66,331 111,736 59,333.20 19,133.11 
   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 r0=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 Q6 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.

Figure 3:

The left panel shows the traces of the BOP values Q4, Q6, and for all minimum configurations found by NM (black lines) and CMA-ES (dots). The trace for the best NM structure is highlighted in gray. The right panel shows this best CK220 configuration of NM and the optimal configuration as found by CMA-ES. The NM structure is dodecahedral with distorted pentagonal faces. One face is highlighted in either configuration.

Figure 3:

The left panel shows the traces of the BOP values Q4, Q6, and for all minimum configurations found by NM (black lines) and CMA-ES (dots). The trace for the best NM structure is highlighted in gray. The right panel shows this best CK220 configuration of NM and the optimal configuration as found by CMA-ES. The NM structure is dodecahedral with distorted pentagonal faces. One face is highlighted in either configuration.

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 CK18 and CK220 cluster problems. For the CK18 clusters, the Q4, Q6, , and values of CMA-ES’ mean converge to the optimal values after about 200 generations. For the CK220 clusters, stable optimal BOP values for Q4, Q6, and are reached after about 2,600 generations. The values (data not shown) do not converge to the optimum CK220 cluster. Such variation has also been observed for other instances. We therefore suggest to use the triplet Q4, Q6, 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).

Figure 4:

Trajectories of a typical CMA-ES run. BOP values of CMA-ES’ mean m(g) are plotted versus generation number (g). The left panel shows a CMA-ES trajectory from the CK18 cluster problem, the right panel one from the CK220 cluster problem. The dots represent the BOP values of the optimal solutions.

Figure 4:

Trajectories of a typical CMA-ES run. BOP values of CMA-ES’ mean m(g) are plotted versus generation number (g). The left panel shows a CMA-ES trajectory from the CK18 cluster problem, the right panel one from the CK220 cluster problem. The dots represent the BOP values of the optimal solutions.

The failure of the Nelder-Mead algorithm on the CK220 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 CK220 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 = 104n. 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 Q4, Q6, and for all observed structures.

We first report the results on the CK1 clusters. For these clusters, assuming an exponential scaling of the energies ECK1 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.

Figure 5:

The left panel shows ECK1 of the putative ground states versus the number of atoms in cluster N. The dots represent the results from the CMA-ES runs, and the dashed curve is the best exponential least-squares fit. The right panel shows box plots of the number of FES needed until CMA-ES converges as a function of N.

Figure 5:

The left panel shows ECK1 of the putative ground states versus the number of atoms in cluster N. The dots represent the results from the CMA-ES runs, and the dashed curve is the best exponential least-squares fit. The right panel shows box plots of the number of FES needed until CMA-ES converges as a function of N.

Table 3:
Summary statistics of the putative CK1 ground states and low-energy minima as identified by CMA-ES. The number of particles, energy, and the BOP values Q4, Q6, and are reported for all instances. For N=10, 14, 16, 22, 23, 25, multiple minima were identified. The putative ground states are labeled “A”, other minima “B” or “C”.
NEnergyQ4Q6
0.1084 −0.0931 
0.4630 0.3750 0.7408 −0.0463 
1.0583 0.5092 0.6285 0.0132 
1.1029 0.5619 0.4369 0.0076 
1.9838 0.6250 0.4556 0.0466 
3.1501 0.7638 0.3536 0.0132 
4.6241 0.5118 0.2861 0.0598 
6.3434 0.5092 0.6285 0.0132 
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.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 
NEnergyQ4Q6
0.1084 −0.0931 
0.4630 0.3750 0.7408 −0.0463 
1.0583 0.5092 0.6285 0.0132 
1.1029 0.5619 0.4369 0.0076 
1.9838 0.6250 0.4556 0.0466 
3.1501 0.7638 0.3536 0.0132 
4.6241 0.5118 0.2861 0.0598 
6.3434 0.5092 0.6285 0.0132 
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.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 CK112 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 Q4=0, Q6=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 CK112 can be found rapidly (in less than FES on average) and robustly (all runs converge to the putative ground state).

For CK114 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).

Figure 6:

Putative ground state configuration (CK114A) and two low-energy stable configurations (CK114B and CK114C) of the CK114 cluster.

Figure 6:

Putative ground state configuration (CK114A) and two low-energy stable configurations (CK114B and CK114C) of the CK114 cluster.

The energy landscape of CK116 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. CK116A consists of two opposite rotated square faces (like the anticube), and triangular faces otherwise. CK116B 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 CK116A) 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 CK116 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.

Figure 7:

Putative ground state configuration (CK116A) and a competing sub-optimal configuration (CK116B) of the CK116 cluster.

Figure 7:

Putative ground state configuration (CK116A) and a competing sub-optimal configuration (CK116B) of the CK116 cluster.

Figure 8:

The left panel shows the -dependent transition probability P for CMA-ES. The right panel depicts a typical trajectory of CMA-ES’ mean in the plane for . Each configuration is color-coded by the .

Figure 8:

The left panel shows the -dependent transition probability P for CMA-ES. The right panel depicts a typical trajectory of CMA-ES’ mean in the plane for . Each configuration is color-coded by the .

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 ECK2 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 CK116 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.

Figure 9:

The left panel shows ECK2 of the putative ground states versus the number of atoms in the cluster N. The dots represent the results from the CMA-ES runs, the dashed curve is the best quadratic least-squares fit. The right panel shows box plots of the number of FES needed until CMA-ES converges as a function of N.

Figure 9:

The left panel shows ECK2 of the putative ground states versus the number of atoms in the cluster N. The dots represent the results from the CMA-ES runs, the dashed curve is the best quadratic least-squares fit. The right panel shows box plots of the number of FES needed until CMA-ES converges as a function of N.

Figure 10:

Putative tetrahedral ground state configuration (CK14A) and a competing suboptimal pyramidal configuration (CK24B) of the CK24 cluster.

Figure 10:

Putative tetrahedral ground state configuration (CK14A) and a competing suboptimal pyramidal configuration (CK24B) of the CK24 cluster.

Table 4:
Summary statistics of the putative CK2 ground states and local minima as identified by CMA-ES. The number of particles, energy, and the BOP values Q4, Q6, and are reported for all instances. For , multiple minima were identified. The putative ground states are labeled “A”.
VEnergyQ4Q6
−0.0931 
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 
6.0977 0.6250 0.4556 0.0466 
12.0076 0.7638 0.3536 0.0132 
29.2253 0.5536 0.0625 −0.0931 
49.3528 0.3736 0.2502 −0.0931 
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.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.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 
VEnergyQ4Q6
−0.0931 
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 
6.0977 0.6250 0.4556 0.0466 
12.0076 0.7638 0.3536 0.0132 
29.2253 0.5536 0.0625 −0.0931 
49.3528 0.3736 0.2502 −0.0931 
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.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.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.

Table 5:
Suggested benchmark settings for the CK1/CK2 cluster test suite.
ProblemsCK1, CK2
Runs/problem 25 
n  
MAX_FES  
Termination If FES=MAX_FES or 
  
Initialization and bounds Uniformly at random in  
ProblemsCK1, CK2
Runs/problem 25 
n  
MAX_FES  
Termination If FES=MAX_FES or 
  
Initialization and bounds Uniformly at random in  

5  Lennard-Jones Clusters

Energy landscapes of collections of atoms that interact according to the Lennard-Jones (LJ) pair potential are among the best studied models in theoretical chemistry and biophysics. In theoretical chemistry, the LJ potential is widely used to model the behavior of noble gases such as argon. Biophysicists use the LJ potential to model hydrophobic forces in polymers such as proteins and alkanes. The problem of finding minimum-energy configurations of LJ clusters fascinated researchers for over three decades and is regularly used as a standard test case for first-order search heuristics. In LJ clusters, each pair of atoms interacts through the following pair potential:
formula
11
where rij=dE(pi, pj) and pi=(xi, yi, zi), the 3D position of atom 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 ELJ of a cluster of N LJ atoms is given by:
formula
12
Again, and rij=dE(pi, pj). 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.
Figure 11:

The Lennard-Jones pair potential uLJ(r) versus distance r for . The minimum is at with energy . For , the potential asymptotically approaches 0.

Figure 11:

The Lennard-Jones pair potential uLJ(r) versus distance r for . The minimum is at with energy . For , the potential asymptotically approaches 0.

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 LJ24 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 LJ38 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=3N−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 xy 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 LJ26 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 Q4, Q6, and with the suggested r0=1.391 for LJ clusters (Doye et al., 1999b). The data are summarized in Table 6. Most of the structures have a low Q4 value, indicating icosahedral symmetry (Steinhardt et al., 1983). As previously mentioned, the LJ13 ground state is identical (with identical BOP values) to the putative CK112 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 LJ19 instance, Kunz and Berry (1995) estimated the number of geometrically distinct local minima to be around half a million.

Table 6:
Summary statistics of the putative LJ ground states from the Cambridge Cluster Database (Wales et al., 2009). The number of particles, energy, and the BOPs Q4, Q6, and are given for each instance.
NEnergyQ4Q6
−3 0.3750 0.7408 −0.0463 
−6 0.1909 0.5745 −0.0132 
−9.1039 0.0013 0.4297 0.0314 
−12.7121 0.1909 0.5745 −0.0132 
−16.5054 0.0148 0.1604 −0.0931 
−19.8215 0.0644 0.1467 0.0015 
−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.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 
NEnergyQ4Q6
−3 0.3750 0.7408 −0.0463 
−6 0.1909 0.5745 −0.0132 
−9.1039 0.0013 0.4297 0.0314 
−12.7121 0.1909 0.5745 −0.0132 
−16.5054 0.0148 0.1604 −0.0931 
−19.8215 0.0644 0.1467 0.0015 
−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.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 
Figure 12:

Putative ground states for LJ clusters with seven atoms (top left), 13 atoms (top right), and 19 atoms (bottom). The packing of LJ13 is a Mackay icosahedron. The putative optimal CK112 cluster, overlaid in black, shares the same structure.

Figure 12:

Putative ground states for LJ clusters with seven atoms (top left), 13 atoms (top right), and 19 atoms (bottom). The packing of LJ13 is a Mackay icosahedron. The putative optimal CK112 cluster, overlaid in black, shares the same structure.

We now provide numerical experiments on three select instances: LJ7, LJ13, and LJ19. 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). LJ13 and LJ19 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 ps for CMA-ES runs on LJ7, LJ13, and LJ19 are summarized in Table 7. CMA-ES always finds the optimal structure of LJ7. In 24 out of the 25 runs, it also solves the LJ13 cluster problem. Inspection of Doye's and Wales's disconnectivity graphs for LJ13 (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 LJ19 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.

Table 7:
Statistics about the number of FES for CMA-ES runs that reach the globally optimal configurations for the LJ7, LJ13, and LJ19 cluster problems, and the success rate ps.
Min7thMedian19thMaxMeanstdps
LJ7 3,805 5,509 15,108 23,980 89,819 18,220.36 18,218.27 
LJ13 13,007 51,827 101,726 178,359 — 109,377.12 80,713.15 0.96 
LJ19 31,726 — — — — 270,870.50 338,201.40 0.08 
Min7thMedian19thMaxMeanstdps
LJ7 3,805 5,509 15,108 23,980 89,819 18,220.36 18,218.27 
LJ13 13,007 51,827 101,726 178,359 — 109,377.12 80,713.15 0.96 
LJ19 31,726 — — — — 270,870.50 338,201.40 0.08 
Table 8:
Suggested benchmark settings for the LJ cluster test suite.
ProblemsLJ
Runs/problem 25 
n  
MAX_FES  
Termination If FES=MAX_FES or 
  
Initialization and bounds Uniformly at random in [−2, 2]n 
ProblemsLJ
Runs/problem 25 
n  
MAX_FES  
Termination If FES=MAX_FES or 
  
Initialization and bounds Uniformly at random in [−2, 2]n 

5.2  The LJ38 Cluster as a High-Dimensional Benchmark with Tunable Landscape Topology

All cluster instances introduced so far have a static landscape topology that is solely determined by the underlying pair potential and the number of atoms. We complement the set of cluster benchmarks with a problem instance that (1) is high-dimensional, but solvable in reasonable time, and (2) exhibits a tunable landscape topology. The benchmark is based on the LJ38 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 LJ38 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 LJ38 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 LJ38 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:
formula
13
where and where pcm 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 Qc(xico)=96.1624; and for the best fcc structure Qc(xfcc)=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 LJ38 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.
Figure 13:

Sketch of the -dependent evolution of the LJ38 energy landscape. The x axis represents a suitable order parameter that can discriminate between different cluster topologies, the y axis represents the potential energy E. The global topology of the landscape is a double funnel with a suboptimal ico structure at the bottom of the wider funnel. The narrow funnel contains the global minimum with fcc symmetry. Increasing gradually changes the landscape topology from a double funnel to a single funnel.

Figure 13:

Sketch of the -dependent evolution of the LJ38 energy landscape. The x axis represents a suitable order parameter that can discriminate between different cluster topologies, the y axis represents the potential energy E. The global topology of the landscape is a double funnel with a suboptimal ico structure at the bottom of the wider funnel. The narrow funnel contains the global minimum with fcc symmetry. Increasing gradually changes the landscape topology from a double funnel to a single funnel.

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 xfcc 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, Q4 is best suited for discriminating between fcc and ico structures, with Q4(xico)=0 and Q4(xfcc)=0.1909. Table 10 summarizes our proposed specification for the tunable LJ38 cluster benchmark. The bounds are far from being tight with respect to the optimal structure. Both xico and xfcc 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.

Table 9:
Characteristics of the putative ground states of compressed LJ38 clusters for different values of . We report the final energy and the BOPs Q4, Q6, and of each structure.
EnergyQ4Q6
−173.92843 0.19090 0.57446 −0.01316 
0.5 −128.42914 0.19090 0.57445 −0.01316 
−83.50595 0.19090 0.57444 −0.01316 
1.5 −39.08437 0.19091 0.57443 −0.01316 
4.89246 0.19091 0.57442 −0.01316 
2.5 48.46946 0.19091 0.57441 −0.01316 
91.68310 0.19092 0.57440 −0.01316 
3.5 134.56362 0.19092 0.57439 −0.01316 
177.13651 0.19092 0.57438 −0.01316 
4.5 219.42358 0.19092 0.57437 −0.01316 
261.44371 0.19092 0.57436 −0.01316 
EnergyQ4Q6
−173.92843 0.19090 0.57446 −0.01316 
0.5 −128.42914 0.19090 0.57445 −0.01316 
−83.50595 0.19090 0.57444 −0.01316 
1.5 −39.08437 0.19091 0.57443 −0.01316 
4.89246 0.19091 0.57442 −0.01316 
2.5 48.46946 0.19091 0.57441 −0.01316 
91.68310 0.19092 0.57440 −0.01316 
3.5 134.56362 0.19092 0.57439 −0.01316 
177.13651 0.19092 0.57438 −0.01316 
4.5 219.42358 0.19092 0.57437 −0.01316 
261.44371 0.19092 0.57436 −0.01316 
Table 10:
Suggested benchmark settings for the tunable LJ38 cluster test case.
ProblemLJ38 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 
ProblemLJ38 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 LJ38 test case without compression. Since many real-world applications involve multi-funnel landscapes, we believe that the tunable LJ38 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 Qc of the LJ38 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

Adip
,
A. B.
(
2005
).
NP-hardness of the cluster minimization problem revisited
.
Journal of Physics A: Mathematical and General
,
38
(
40
):
8487
,
doi: 10.1088/0305-4470/38/40/001
Auger
,
A.
, and
Hansen
,
N
. (
2005
).
A restart CMA evolution strategy with increasing population size
. In
Proceedings of the IEEE Congress on Evolutionary Computation (CEC 2005)
, Vol.
2
, pp.
1769
1776
.
Barron
,
C.
,
Gómez
,
S.
, and
Romero
,
D
. (
1996
).
Archimedean polyhedron structure yields a lower energy atomic cluster
.
Applied Mathematics Letters
,
9
(
5
):
75
78
.
Born
,
M.
, and
Oppenheimer
,
R
. (
1927
).
Zur quantentheorie der molekeln
.
Annalen der Physik
,
389
(
20
):
457
484
.
Call
,
S. T.
,
Zubarev
,
D. Y.
, and
Boldyrev
,
A. I
. (
2007
).
Global minimum structure searches via particle swarm optimization
.
Journal of Computational Chemistry
,
28
(
7
):
1177
1186
.
Catlow
,
C. R. A.
,
Bromley
,
S. T.
,
Hamad
,
S.
,
Mora-Fonz
,
M.
,
Sokol
,
A. A.
, and
Woodley
,
S. M
. (
2010
).
Modelling nano-clusters and nucleation
.
Physical Chemistry Chemical Physics
,
12
(
4
):
786
811
.
Cleri
,
F.
, and
Rosato
,
V
. (
1993
).
Tight-binding potentials for transition metals and alloys
.
Physical Review B
,
48
(
1
):
22
33
.
Cohn
,
H.
, and
Kumar
,
A
. (
2009
).
Algorithmic design of self-assembling structures
.
Proceedings of the National Academy of Sciences of the United States of America
,
106
(
24
):
9570
9575
.
Cox
,
G.
,
Berry
,
R. S.
, and
Johnston
,
R. L
. (
2006
).
Characterizing potential surface topographies through the distribution of saddles and minima
.
Journal of Physical Chemistry A
,
110
(
40
):
11543
11550
.
Deaven
,
D.
,
Tit
,
N.
,
Morris
,
J.
, and
Ho
,
K.
(
1996
).
Structural optimization of Lennard-Jones clusters by a genetic algorithm
.
Chemical Physical Letters
,
256
(
1–2
):
195
200
.
Doye
,
J.
(
2000
).
Effect of compression on the global optimization of atomic clusters
.
Physical Review E
,
62
(
6, Part B
):
8753
8761
.
Doye
,
J. P. K.
,
Miller
,
M. A.
, and
Wales
,
D. J
. (
1999a
).
Evolution of the potential energy surface with size for Lennard-Jones clusters
.
Journal of Chemical Physics
,
111
(
18
):
8417
8428
.
Doye
,
J. P. K.
,
Miller
,
M. A.
, and
Wales
,
D. J
. (
1999b
).
The double-funnel energy landscape of the 38-atom Lennard-Jones cluster
.
Journal of Chemical Physics
,
110
(
14
):
6896
6906
.
Finck
,
S.
,
Hansen
,
N.
,
Ros
,
R.
, and
Auger
,
A.
(
2009
).
Real-parameter black-box optimization benchmarking 2009: Presentation of the noiseless functions contents
, http://citesecrx.1st.psu.edu-viewdoc-summary?doi=10.1.140.4404
Finnis
,
M. W.
, and
Sinclair
,
J. E
. (
1984
).
A simple empirical N-body potential for transition metals
.
Philosophical Magazine A
,
50
(
1
):
45
55
.
Hansen
,
N.
(
2008
).
The CMA evolution strategy: A tutorial
.
Retrieved from
http://www.lri.fr/~hansen/cmatutorial.pdf
Hansen
,
N
. (
2009a
).
Benchmarking a bi-population CMA-ES on the BBOB-2009 noisy testbed
. In
Proceedings of the 11th Annual Conference Companion on Genetic and Evolutionary Computation Conference: Late Breaking Papers, GECCO ’09
, pp.
2397
2402
.
Hansen
,
N.
(
2009b
).
Benchmarking the Nelder-Mead downhill simplex algorithm with many local restarts
. In
Proceedings of the 11th Annual Conference Companion on Genetic and Evolutionary Computation Conference
,
GECCO ’09
, pp.
2403
2408
.
Hartke
,
B.
(
2004
).
Application of evolutionary algorithms to global cluster geometry optimization
. In
R. L. Johnston (Ed.)
,
Applications of Evolutionary Computation in Chemistry
, Vol.
110
(pp.
33
53
).
Berlin
:
Springer
.
Hasse
,
R. W
. (
1991
).
Structure and magic numbers of large Lennard-Jones quasicrystals and crystals
.
Physics Letters A
,
161
(
2
):
130
134
.
Hoare
,
M. R
. (
1979
).
Structure and dynamics of simple microclusters
.
New York
:
Wiley
.
Kunz
,
R.
, and
Berry
,
R
. (
1995
).
Statistical interpretation of topographies and dynamics of multidimensional potentials
.
Journal of Chemical Physics
,
103
(
5
):
1904
1912
.
Leary
,
R
. (
2000
).
Global optimization on funneling landscapes
.
Journal of Global Optimization
,
18
(
4
):
367
383
.
Lennard-Jones
,
J
. (
1924
).
On the determination of molecular fields, II, From the equation of state of a gas
.
Proceedings of the Royal Society of London Series A, Containing Papers of a Mathematical and Physical Character
,
106
(
738
):
463
477
.
Li
,
Z.
, and
Scheraga
,
H. A
. (
1987
).
Monte Carlo-minimization approach to the multiple-minima problem in protein folding
.
Proceedings of the National Academy of Sciences of the United States of America
,
84
(
19
):
6611
6615
.
Meng
,
G.
,
Arkus
,
N.
,
Brenner
,
M. P.
, and
Manoharan
,
V. N
. (
2010
).
The free-energy landscape of clusters of attractive hard spheres
.
Science
,
327
(
5965
):
560
563
.
Morse
,
P
. (
1929
).
Diatomic molecules according to the wave mechanics, II; Vibrational levels
.
Physical Review
,
34
(
1
):
57
64
.
Müller
,
C. L.
,
Baumgartner
,
B.
, and
Sbalzarini
,
I. F
. (
2009
).
Particle swarm CMA evolution strategy for the optimization of multi-funnel landscapes
. In
Proceedings of the IEEE Congress on Evolutionary Computation (CEC 2009)
, pp.
2685
2692
.
Müller
,
C. L.
, and
Sbalzarini
,
I. F.
(
2009
).
A tunable real-world multi-funnel benchmark problem for evolutionary optimization—And why parallel island models might remedy the failure of CMA-ES on it
. In
International Joint Conference on Computational Intelligence, IJCCI
, pp.
248
253
.
Müller
,
S. D.
,
Schraudolph
,
N. N.
, and
Koumoutsakos
,
P
. (
2003
).
Evolutionary and gradient-based algorithms for Lennard-Jones cluster optimization
. In
Proceedings of the 5th Annual Conference on Genetic and Evolutionary Computation, GECCO ’03
.
Nelder
,
J. A.
, and
Mead
,
R
. (
1965
).
A simplex method for function minimization
.
The Computer Journal
,
7
(
4
):
308
313
.
Northby
,
J
. (
1987
).
Structure and binding of Lennard-Jones clusters: 13 147
.
Journal of Chemical Physics
,
87
(
10
):
6166
6177
.
Osada
,
R.
,
Funkhouser
,
T.
,
Chazelle
,
B.
, and
Dobkin
,
D
. (
2002
).
Shape distributions
.
ACM Transactions on Graphics
,
21
(
4
):
807
832
.
Rechtsman
,
M. C.
,
Stillinger
,
F. H.
, and
Torquato
,
S
. (
2005
).
Optimized interactions for targeted self-assembly: Application to a honeycomb lattice
.
Physical Review Letters
,
95
(
22
):
228301
.
Rechtsman
,
M. C.
,
Stillinger
,
F. H.
, and
Torquato
,
S
. (
2006a
).
Designed interaction potentials via inverse methods for self-assembly
.
Physical Review E
,
73
(
1
):
011406
.
Rechtsman
,
M. C.
,
Stillinger
,
F. H.
, and
Torquato
,
S
. (
2006b
).
Self-assembly of the simple cubic lattice with an isotropic potential
.
Physical Review E
,
74
(
2
):
021404
.
Rechtsman
,
M. C.
,
Stillinger
,
F. H.
, and
Torquato
,
S
. (
2007
).
Synthetic diamond and wurtzite structures self-assemble with isotropic pair interactions
.
Physical Review E
,
75
(
3
):
031403
.
Sloane
,
N.
,
Hardin
,
R.
,
Duff
,
T.
, and
Conway
,
J
. (
1995
).
Minimal-energy clusters of hard spheres
.
Discrete and Computational Geometry
,
14
(
1
):
237
259
.
Sloane
,
N. J. A.
,
Hardin
,
R. H.
,
Duff
,
T. S.
, and
Conway
,
J. H.
(
1997
).
The sphere packing cluster database
.
Information Sciences Research, AT&T Shannon Lab, Florham Park, New Jersey
.
Retrieved from
http://www2.research.att.com/~njas/cluster/index.html
Sloane
,
N. J. A.
,
Hardin
,
R. H.
, and
Smith
,
W. D.
(
2000
).
Spherical codes: Nice arrangements of points on a sphere in various dimensions. Information Sciences Research, AT&T Shannon Lab, Florham Park, New Jersey
.
Retrieved from
http://www2.research.att.com/~njas/packings/
Smit
,
S.
, and
Eiben
,
A.
(
2010
).
Beating the “world champion” evolutionary algorithm via REVAC tuning
. In
Proceedings of the 2010 IEEE Congress on Evolutionary Computation (CEC)
, pp.
1
8
.
Steinhardt
,
P. J.
,
Nelson
,
D. R.
, and
Ronchetti
,
M
. (
1983
).
Bond-orientational order in liquids and glasses
.
Physical Review B
,
28
(
2
):
784
805
.
Stillinger
,
F. H.
, and
Weber
,
T. A
. (
1984
).
Packing structures and transitions in liquids and solids
.
Science
,
225
(
4666
):
983
989
.
Suganthan
,
P. N.
,
Hansen
,
N.
,
Liang
,
J. J.
,
Deb
,
K.
,
Chen
,
Y.-P.
,
Auger
,
A.
, and
Tiwari
,
S.
(
2005
).
Problem definitions and evaluation criteria for the CEC 2005 special session on real-parameter optimization
.
Technical report, Nanyang Technological University, Singapore (Tech. Rep. 2009005)
.
Sutton
,
A. P.
, and
Chen
,
J
. (
1990
).
Long-range Finnis-Sinclair potentials
.
Philosophical Magazine Letters
,
61
(
3
):
139
146
.
Thomson
,
J. J
. (
1904
).
On the structure of the atom: An investigation of the stability and periods of oscillation of a number of corpuscles arranged at equal intervals around the circumference of a circle; With application of the results to the theory of atomic structure
.
Philosophical Magazine
,
7
(
39
):
237
265
.
Torquato
,
S
. (
2009
).
Inverse optimization techniques for targeted self-assembly
.
Soft Matter
,
5
(
6
):
1157
1173
.
Van Hoyweghen
,
C.
, and
Naudts
,
B
. (
2000
).
Symmetry in the search space
. In
Proceedings of the IEEE Congress on Evolutionary Computation (CEC 2000)
, Vol.
2
, pp.
1072
1078
.
Van Hoyweghen
,
C.
,
Naudts
,
B.
, and
Goldberg
,
D
. (
2002
).
Spin-flip symmetry and synchronization
.
Evolutionary Computation
,
10
(
4
):
317
344
.
Wales
,
D
. (
2004
).
Energy landscapes: Applications to clusters, biomolecules and glasses
.
Cambridge, UK
:
Cambridge University Press
.
Wales
,
D. J
. (
2005
).
Energy landscapes and properties of biomolecules
.
Physical Biology
,
2
(
4
):
S86
S93
.
Wales
,
D. J.
, and
Doye
,
J. P. K
. (
1997
).
Global optimization by basin-hopping and the lowest energy structures of Lennard-Jones clusters containing up to 110 atoms
.
Journal of Physical Chemistry A
,
101
(
28
):
5111
5116
.
Wales
,
D. J.
,
Doye
,
J. P. K.
,
Dullweber
,
A.
,
Hodges
,
M. P.
,
Naumkin
,
F. Y.
,
Calvo
,
F.
,
Hernández-Rojas
,
J.
, and
Middleton
,
T. F.
(
2009
).
The Cambridge cluster database. Retrieved from
http://www-wales.ch.cam.ac.uk/CCD.html
Wales
,
D. J.
, and
Scheraga
,
H.
(
1999
).
Review: Chemistry—Global optimization of clusters, crystals, and biomolecules
.
Science
,
285
(
5432
):
1368
1372
.
Weber
,
T. A.
, and
Stillinger
,
F. H
. (
1985
).
Interactions, local order, and atomic-rearrangement kinetics in amorphous nickel-phosphorous alloys
.
Physical Review B
,
32
(
8
):
5402
5411
.
Weisstein
,
E. W.
(
2010
).
Wigner 3j-symbol. From MathWorld—A Wolfram Web Resource
.
Retrieved from
http://mathworld.wolfram.com/Wigner3j-Symbol.html
Whitley
,
D.
(
2010
).
Personal communication. Chair, Department of Computer Science, Colorado State University, Fort Collins
.
Wille
,
L
. (
1987
).
Minimum-energy configurations of atomic clusters: New results obtained by simulated annealing
.
Chemical Physics Letters
,
133
(
5
):
405
410
.
Wille
,
L. T.
, and
Vennik
,
J
. (
1985
).
Computational complexity of the ground-state determination of atomic clusters
.
Journal of Physics A: Mathematical and General
,
18
(
8
):
L419
.
Wright
,
S.
(
1932
).
The roles of mutation, inbreeding, crossbreeding, and selection in evolution
.
Proceedings of the Sixth International Congress on Genetics
, pp.
356
366
.
Wu
,
C. W.
, and
Verma
,
D
. (
2008
).
A sensor placement algorithm for redundant covering based on Riesz energy minimization
. In
International Symposium on Circuits and Systems, ISCAS
, pp.
2074
2077
.

Note