Global Optimization by Basin-Hopping and the Lowest Energy Structures of Lennard-Jones

Global Optimization by Basin-Hopping and the Lowest Energy Structures of...

(Parte 1 de 3)

Global Optimization by Basin-Hopping and the Lowest Energy Structures of Lennard-Jones Clusters Containing up to 110 Atoms

David J. Wales* UniVersity Chemical Laboratories, Lensfield Road, Cambridge CB2 1EW, U.K.

Jonathan P. K. Doye FOM Institute for Atomic and Molecular Physics, Kruislaan 407, 1098 SJ Amsterdam, The Netherlands

ReceiVed: March 19, 1997; In Final Form: April 29, 1997X

We describe a global optimization technique using “basin-hopping” in which the potential energy surface is transformed into a collection of interpenetrating staircases. This method has been designed to exploit the features that recent work suggests must be present in an energy landscape for efficient relaxation to the global minimum. The transformation associates any point in configuration space with the local minimum obtained by a geometry optimization started from that point, effectively removing transition state regions from the problem. However, unlike other methods based upon hypersurface deformation, this transformation does not change the global minimum. The lowest known structures are located for all Lennard-Jones clusters up to 110 atoms, including a number that have never been found before in unbiased searches.

I. Introduction

Global optimization is a subject of intense current interest.1

Improved global optimization methods could be of great economic importance, since improved solutions to traveling salesman-type problems, the routing of circuitry in a chip, the active structure of a biomolecule, etc., equate to reduced costs or improved performance. In chemical physics the interest in efficient global optimization methods stems from the common problemof findingthe lowestenergyconfigurationof a (macro)- molecular system. For example, it seems likely that the native structure of a protein is structurally related to the global minimum of its potential energy surface (PES). If this global minimum could be found reliably from the primary amino acid sequence, this knowledge would provide new insight into the nature of protein folding and save biochemists many hours in the laboratory. Unfortunately, this goal is far from being realized. Instead the development of global optimization methods has usually concentrated on much simpler systems.

Lennard-Jones (LJ) clusters represent one such test system. Here the potential is

where and 21/6ó are the pair equilibrium well depth and separation, respectively. We will employ reduced units, i.e., ) ó ) 1 throughout. Much of the initial interest in LJ clusters was motivated by a desire to calculate nucleation rates for noble gases. However, as a result of the wealth of data generated, the LJ potential has been used not only for studying global optimization but also the effects of finite size on phase transitions such as melting. Through the combined efforts of many workers, likely candidates for the global minima of LJN clusters have been found up to N ) 147.2-16 This represents a significantachievementsince extrapolationof Tsai and Jordan’s comprehensive enumeration of minima for small LJ clusters17 suggests that the PES of the 147-atom cluster possesses of the order of 1060 minima.18

at N ) 13, 5, 147,At most intermediate sizes the global

Previous studies have revealed that the Mackay icosahedron19 provides the dominant structural motif for LJ clusters in the size range of 10-150 atoms. Complete icosahedra are possible minimum consists of a Mackay icosahedron at the core covered by a low-energy overlayer. As a consequence of the phase behavior of LJ clusters, finding these global minima is relatively easy. Studies have shown that in the region of the solid-liquid transition the cluster is observed to change back and forth between a liquid-like form and icosahedral structures.20 As a result of this “dynamic coexistence,” a method as crude as molecular dynamics within the melting region coupled with systematic minimization of configurations generated by the trajectory is often sufficient to locate the global minimum.21

However, there are a number of sizes at which the global minimum is not based on an icosahedral structure. These clusters are illustrated in Figure 1. For LJ38 the lowest energy structure is a face-centered-cubic(fcc) truncated octahedron,13,14 and for N ) 75, 76, 7, 102, 103, and 104, geometries based on Marks’decahedra22 are lowestin energy.14,15 For these cases, finding the lowest minimum is much harder because the global minimumof freeenergyonlybecomesassociatedwiththe global potential energy minimum at temperatures well below melting where the dynamics of structural relaxation are very slow. For

LJ38, the microcanonicaltemperaturefor the transitionfrom facecentered cubic to icosahedral structures has been estimated to be about 0.12 k-1, where k is the Boltzmann constant, and forX Abstract published in AdVance ACS Abstracts, June 15, 1997.

Figure 1. Nonicosahedral Lennard-Jones global minima.

S1089-5639(97)0984-5 C: $14.0 © 1997 American Chemical Society

LJ75 the estimate for the decahedral to icosahedral transition is about 0.09 k-1.23 (For comparison, melting typically occurs

The topography of the PES can also play a key role in determining the ease of global optimization.24 A detailed study of the LJ38 PES has shown that there is a large energy barrier between the fcc and icosahedral structures,25 which correspond to well-separated regions of the PES. Furthermore, fcc and decahedral structures have less polytetrahedral character than icosahedralstructures,and hence they have less in common with the liquid-like state, which is characterized by disordered polytetrahedral packing.26,27 Since the vast majority of configuration space is dominated by “liquid-like” configurations, it is therefore harder to find global minima based upon fcc and decahedral packing using unbiased searches.

These considerations explain why global optimization methods have only recently begun to find the truncated octahedron13,16,28,29 and why, until now, the Marks’ decahedron has never been found by an unbiased global optimization method.

The greater difficulty of finding the LJ75 global minimum compared to LJ38 can probably be explained by the slightly smaller transition temperature, the sharper transition,23 and the much larger number of minima on the LJ75 PES.

Before we consider the effectiveness of different global optimizationmethodsfor Lennard-Jonesclusters,it is interesting to note that the use of physical principles to construct good candidate structures for the global minima2-5,8,14,15 or to reduce the configuration space that needs to be searched10-12 led to the initial discovery of 93% of the LJ global minima in this size range. It seems that physicalinsight into a specificproblem will often be able to beat unbiased global optimization, a view expressed by Ngo et al.30 in their discussion of computational complexity.

One difficulty in evaluating the relative performances of different global optimization methods is that, too often, the methods have only been applied to small clusters, or to larger clusters with global minima that are especially stable, such as

LJ55. It is also difficult to draw any firm conclusions about how efficient different methods may be when the number of searchesemployedvarieswidely. However,it seems reasonable to suggest two hurdles that any putative global optimization approach should aspire to. The first is the location of the truncated octahedron for LJ38, and any method which fails this test is unlikely to be useful. The second hurdle is the location of the Marks’ decahedron for LJ75; this problem poses a much more severe test for an unbiased search and one which does not appear to have been passed until the present work.

The most successfulglobaloptimizationresultsfor LJ clusters reported to date are for genetic algorithms.16,28,31,32 These methods mimic some aspects of biological evolution: a population of clusters evolves to low energy by mutation and mating of structures, along with selection of those with low potential energy. To be successful, new configurations produced by “genetic manipulation” are mapped onto minima by a local optimization algorithm such as the conjugate gradient method. The study by Deaven et al. is particularlyimpressive,since these workers matched or beat all the lowest energy minima that they knew of up to N ) 100, including the truncated octahedron (although they probably missed the global minima at N ) 69 and 75-78). Niesse and Mayne were also able to locate the

LJ38 truncated octahedron, and report that this structure took about 25 times longer to find than the icosahedralglobal minima of the neighboring sizes.

Another class of global optimization techniques, sometimes called hypersurface deformation methods, attempts to simplify the problem by applying a transformation to the PES which smoothes it and reduces the number of local minima.3,34 The global minimum of the deformed PES is then mapped back to the original surface in the hope that this will lead back to the global minimum of the original PES. The distinctions between the variousmethodsof this type lie in the type of transformations that are used, which include applying the diffusion equation,35 increasing the range of the potential13,36 and shifting the position of the potential minimum toward the origin.37 The performance of hypersurfacedeformationmethodshas been variable: Pillardy and Piela13 managed to find the 38-atom truncated octahedron, but other workers report difficulties35 for the trivial cases of

LJ8 and LJ9 where there are only 8 and 21 minima on the PES, respectively.

Althoughintuitivelyappealing,the problemwith hypersurface deformation is that there is no guarantee that the global minimum on the deformed PES will map onto the global minimum of the original surface. This difficulty is clearly illustrated when we consider Stillinger and Stillinger’s suggestion of increasing the range of the potential:36 it has been shown that the global minimum may in fact depend rather sensitively on the range of the potential, with the appearance of numerous “range-induced” transitions.14,38

Other methods include those based on “annealing”. Such approaches take advantage of the simplification in the free energy landscape that occurs at high temperatures and attempt to follow the free energy global minimum as the temperature is decreased. At 0 K, the free energy global minimum and the global minimum of the PES must coincide. Standard simulated annealing39 was used by Wille to find a few new minima at small sizes9 but does not appear to have been systematically applied to LJ clusters. More sophisticated variants of this techniqueincludegaussiandensityannealingand analogues,40-43 but again some appear to fail at small sizes.42,43

The difficulty with the annealing approach methods is that, if the free energy global minimum changes at low temperatures where dynamical relaxation is slow, the algorithms will become stuck in the structure correspondingto the high temperature free energy global minimum. Such methods are therefore likely to experience difficulties in finding the global minima for LJ38 and LJ75. In the language employed in recent protein-folding literature,4 annealing will fail when Tf < Tg, where Tf is the “folding” temperature below which the global potential energy and free energy minima coincide, and Tg is the “glass” temperature at which the system effectively becomes trapped in a local minimum.

Another method which attempts to reduce the effects of barriers on the PES makes use of quantum tunneling. The diffusion Monte Carlo approach is used to find the ground state wave function, which should become localized at the global minimumas p is decreasedto zero.45 A more rigorousapproach has been applied by Maranas and Floudas, who found upper and lower bounds for the energy of the global minimum. However, the computational expense of this method, which scales as 2N with the number of atoms, means that it has only been used for small systems.46,47 Most of the above studies, along with the recently described “pivot method”48-50 and “taboo search”,51,52 have yet to prove their usefulnessby passing the first hurdle for LJ38 suggested above. However, this does not necessarilymean that theseapproachesshouldbe discounted, since some authors have only applied their algorithmsto smaller clusters and may not have run enough searches to achieve convergence.

In the present work we present the results of a “basinhopping” global optimization technique for Lennard-Jones

5112 J. Phys. Chem. A, Vol. 101, No. 28, 1997 Wales and Doye

clusters. All the known lowest energy structures up to N ) 110 have been located successfully, including three minima not previously reported. (See Tables 1 and 2.) This method is also the first unbiased algorithm to find the global minima based on the Marks decahedron around LJ75 and LJ102. For reference, we collect the rather scattered results previously reported for

LJ clusters to provide a complete catalog of the energies and point groups of the lowest energy minima that we know of. The results have been collected in the first entry of the Cambridge Cluster Database at

I. Method

The present approach has been guided by previous work on energy landscapes which has identified features that enable the system to locate its global minimum efficiently.4 In particular, analysis of model energy landscapes, using a master equation approach for the dynamics, has provided good evidence that such a surface should have a large potential energy gradient and the lowestpossibletransitionstate energiesor rearrangement barriers.53 These results immediately suggest a simple way to transform the PES which does not change the global minimum, nor the relative energies of any local minima. We consider the transformed energy E÷ defined by where X represents the 3N-dimensional vector of nuclear coordinates and min signifies that an energy minimization is performed starting from X. In the present work, energy

TABLE 1: Global Minima of LJN for N e 110 N point group energy/ refa N point group energy/ refa a The reference in which each minimum was first reported (to the best of our knowledge) is given. We intend to maintain an updated database of energies and coordinates for LJ and Morse clusters on our web site:

Lowest Energy Structures of Lennard-Jones Clusters J. Phys. Chem. A, Vol. 101, No. 28, 1997 5113

minimizations were performed using the Polak-Ribiere variant of the conjugate gradient algorithm.54 Hence the energy at any point in configuration space is assigned to that of the local minimum obtained by the given geometry optimization technique, and the PES is mapped onto a set of interpenetrating staircases with plateaus corresponding to the set of configurations which lead to a given minimum after optimization. A schematic view of the staircase topography that results from this transformation is given in Figure 2. These plateaus, or basins of attraction, have been visualized in previous work as a means to compare the efficiency of different transition state searching techniques.5,56

The energy landscape for the function E÷(X) was explored using a canonical Monte Carlo simulation at a constant reduced temperatureof 0.8. At each step, all coordinateswere displaced by a random number in the range [-1,1] times the step size, which was adjusted to give an acceptance ratio of 0.5. The nature of the transformed surface allowed relatively large step sizes of between 0.36-0.40. For each cluster in the range considered, seven separate runs were conducted. Five of these each consistedof 5000 Monte Carlo steps startingfrom different randomlygeneratedconfigurationsof atomsconfinedto a sphere of radius 5.5 reduced units. The subsequentgeometryoptimizations employed a container of radius one plus the value required to contain the same volume per atom as the fcc primitive cell. The container should have little effect on any of our results and is only required to prevent dissociationduring the conjugate gradient optimizations.

The convergence criterion employed for the conjugate gradient optimizations used in the Monte Carlo moves need not be very tight. In the present work we required the root-meansquare (RMS) gradient to fall below 0.01 in reduced units and the energy to change by less than 0.1 between consecutive steps in the conjugate gradient search. Initially it appeared that a tolerance of 0.1 for the RMS gradient was satisfactory, but this was subsequently found to cause problems for clusters containing more than about 60 atoms. The lowest energy structures obtained during the canonical simulation were saved and reoptimized with tolerances of 10-4 and 10-9 for the RMS force and the energy difference,respectively. The final energies are accurate to about six decimal places.

Several other techniqueswere employedin these calculations, namely seeding, freezing and angular moves. Here we used the pair energy per atom, E(i), defined as

so that the total energy is

If the highest pair energy rose above a fraction R of the lowest pair energy then an angular move was employed for the atom in question with all other atoms fixed. R was adjusted to give an acceptance ratio for angular moves of 0.5 and generally convergedto between0.40 and 0.4. Each angulardisplacement consisted of choosing random ı and spherical polar coordinates for the atom in question, taking the origin at the center of mass and replacing the radius with the maximum value in the cluster.

The two remaining runs for each size consisted of only 200

Monte Carlo steps starting from the global minima obtained for the clusters containing one more and one less atom. When starting from LJN-1 the N - 1 atoms were frozen for the first 100 steps, during which only angular moves were attempted for the remaining atom, starting from a random position outside the core. When starting from LJN+1 the atom with the highest pair energy E(i) was removed and 200 unrestrictedMonte Carlo moves were attempted from the resulting geometry.

The above basin-hopping algorithm shares a common philosophy with our previous approach in which steps were taken directly between minima using eigenvector-following to calculate pathways.25 The latter method is similar to that described recently by Barkema and Mousseau57 in their search for wellrelaxed configurations in glasses. Although the computational expense of transition state searches probably makes this method uncompetitive for global optimization, our study illustrated the possible advantage of working in a space in which only the minima are present. The basin-hoppingalgorithmdiffers in that it is applied in configuration space to a transformed surface, rather than in a discrete space of minima, and steps are taken stochastically. The genetic algorithms described by Deaven et al.16 and Niesse and Mayne28 used conjugategradientminimization to refine the local minima which comprise the population of structures that are evolved in their procedure. Hence these authors are in effect studying the same transformed surface as described above, but explore it in a rather different manner. We suspect that the success of their methods is at least partly due to the implicit use of the transformed surface E÷.

The present approach is basically the same as the “Monte

Carlo-minimization” algorithm of Li and Scheraga,58 who applied it to search the conformationalspace of the pentapeptide [Met5]enkephalin. A similar method has recently been used by Baysal and Meirovitch59 to search the conformational space of cyclic polypeptides.

I. Results

(Parte 1 de 3)