The effect of compression on the global optimization of atomic clusters

9 Pages·1.5 MB·English

by  Jonathan Doye

The effect of compression on the global optimization of atomic clusters
Jonathan P. K. Doye
University Chemical Laboratory, Lensfield Road, Cambridge CB2 1EW, UK
(February 1, 2008)
Recently, Locatelli and Schoen proposed a transformation of the potential energy that aids the
global optimization of Lennard-Jonesclusters with non-icosahedral global minima. These cases are
particularlydifficulttooptimizebecausethepotentialenergysurfacehasadoublefunneltopography
withtheglobalminimumatthebottomofthenarrowerfunnel. Hereweanalysetheeffectofthistype
of transformation on the topography of the potential energy surface. The transformation, which
physically corresponds to a compression of the cluster, firstly reduces the number of stationary
points on the potential energy surface. Secondly, we show that for a 38-atom cluster with a face-
centred-cubic global minimum the transformation causes the potential energy surface to become
0
increasingly dominated bythefunnelassociated with theglobal minimum. Thetransformation has
0
been incorporated in thebasin-hoppingalgorithm using a two-phase approach.
0
2
I. INTRODUCTION particularly easy.
n
However,therearesomesizesforwhichtheglobalmin-
a
J imumisnoticosahedral. AtN=38theglobalminimumis
One of the most important types of global optimiza- aface-centred-cubic(fcc)truncatedoctahedron8–10 (38A
6 tion problem, and one which is particularly of interest
in Figure 1), at N=75–77 and 102–104 the global min-
1 to chemical physicists, is the determination of the low- ima are based on Marks26 decahedra10,11 (e.g. 75A in
estenergyconfigurationofamolecularsystem,suchasa
v Figure 1), and at N=98 the global minimum is a Leary
protein,a crystalor acluster.1 However,sucha task can
6 tetrahedron17 (98AinFigure1). ForthesesizesthePES
6 be very difficult because of the large number of minima
has a fundamentally different character. As well as the
0 thatapotentialenergysurface(PES)canhave—itisgen-
wide funnel leading down to the low-energy icosahedral
1 erally expected that the number of minima of a system
structures, there is a much narrower funnel which leads
0 willincreaseexponentiallywithsize.2 Therefore,ifappli-
0 down to the global minimum.22,27 Relaxation down the
cationsto largesystems with realisticdescriptions ofthe
0 PESismuchmorelikelytotakethesystemintothewider
interatomicinteractionsaretobe feasible,itisnecessary
/ funnel where it is then trapped. The time scale for in-
t thatefficientglobaloptimizationalgorithms,whichscale
a terfunnelequilibrationisveryslow28 becauseofthelarge
m well with system size, are developed. energy22 andfreeenergy27 barriersbetweenthe twofun-
A key partofthis developmentis understanding when
- nels.
d and why an algorithm is likely to succeed or fail, be-
n cause, as well as providing useful information about the
o limitations of an algorithm, this physical insight might
c 34A 38A 75A 98A
be utilisedinthedesignofbetter algorithms. This isthe
:
v motivation behind the current paper. Here, we analyse
i the reasons for the success of a recent algorithm when
X
appliedtotheglobaloptimizationofLennard-Jones(LJ)
r
a clusters for some particularly difficult sizes.
34H 38C 75C 98B
The global optimization of LJ clusters has probably
becomethemostcommonbenchmarkforconfigurational
optimization problems.1,3 Putative global minima have
been obtained for all sizes up to 309 atoms,4–17 and up-
to-date databases of these structures are maintained on
the web.18,19 There are two types of difficulty for the LJ
FIG.1. Theglobalminimaandsomelow-lyingminimaof
cluster problem. First, there is the general increase in
the number of minima with cluster size.20,21 Second, on LJ34,LJ38,LJ75andLJ98. 34A,38C,75Cand98Barebased
on Mackay icosahedra. 34H and 98A are Leary tetrahedra.
top of this effect there are size-specific effects related to
38A is a face-centred-cubic truncated octahedron, and 75A
the topography of the PES.22
is a Marks decahedron. The letter gives the energetic rank
For most of the clusters the topography of the PES
oftheminimum,i.e. globalminimaarelabelledwithan‘A’,
aids global optimization. There is a funnel23,24 from the
etc.
high-energyliquid-likeclusterstothelowenergyminima
with structures based up on the Mackay25 icosahedra. As a result these eight clusters are hard to optimize,
When there is a dominant low-energy icosahedral mini- the larger examples being virtually impossible to opti-
mumatthebottomofthefunnel,suchaswhencomplete mize by traditional approaches, such as simulated an-
Mackayicosahedracanbe formed,globaloptimizationis nealing. However, these cases are solvable by a set
1
of methods in which the ‘basin-hopping’ transformation term proportional to r which penalizes long pair
i<j ij
is applied to the PES.13 This transformation is used distances.41 Here, we use a slightly different form, which
by the Monte Carlo minimization29 or basin-hopping again acts to compresPs the cluster. The energy for such
algorithm,13 andimplicitly by allthe mostsuccessfulge- a compressed Lennard-Jones (CLJ) cluster is given by
netic algorithms.12,30–36 The transformation of the PES
works by changing the thermodynamics of the clusters |ri−rc.o.m.|2
E =E + µ , (2)
suchthatthesystemisnowabletopassbetweenthefun- CLJ LJ comp σ2
nelsmoreeasily.37,38However,thenon-icosahedralglobal Xi
minima still take much longer to find than the icosahe- where µcomp is a parameter that determines the magni-
dral global minima,1 and there is no way of knowing if tude of the compressionacting on the cluster, and rc.o.m
one has waited long enough to rule out the possibility is the position of the centre of mass of the cluster. We
of a non-icosahedralglobal minimum. This is illustrated found the additional term to be approximately propor-
by the Leary tetrahedron at N=98. Despite the fact tional to Schoen and Locatelli’s expression, and so the
that powerful optimization techniques had been applied effectofthe twotransformationsonthe PEStopography
to LJ ,13,35,36 the global minimum was discovered only are virtually identical.
98
very recently.17 Subsequently, it was confirmed that this To map the PES topography of these CLJ clusters we
minimum could be found by some of the previously ap- usethesamemethodsasthosewehaveappliedtoLJ22,27
plied methods.39,40 andMorse43clusterstoobtainlargesamplesofconnected
Given this background, it would be useful to develop minima and transition states that provide good repre-
techniquesthataremoreefficientforthesedouble-funnel sentations of the low-energy regions of the PES. The
examples. Two potential approaches have very recently approach involves repeated applications of eigenvector-
been put forward. First, Hartke has achieved improve- following44 to find new transition states and the minima
ments in the genetic algorithm approach by forcing the they connect.
system to maintain a diversity of structural types in In the basin-hopping algorithm,13,45 the transformed
thepopulation,thuspreventingthepopulationbecoming potential energy is given by
concentratedintheicosahedralfunnel.36 Second,Schoen E˜(x)=min{E(x)}, (3)
and Locatelli noted that the exceptions to the icosahe-
dralstructuralmotifsareusuallymoresphericalthanthe where x representsthe vector ofnuclearcoordinatesand
competing icosahedral structures. This is because the min signifies that an energy minimization is performed
exceptionsgenerallyoccuratsizeswhereboth aparticu- starting from x. Hence the energy at any point in con-
larly stable formfor the alternativemorphologyis possi- figurationspaceisassignedto thatofthe localminimum
ble and the icosahedralstructures involve an incomplete obtainedby the minimization, and the transformedPES
overlayer. Therefore,Schoen and Locatelli added a term consists of a set of plateaus or steps each corresponding
to the potential energy favouring compact clusters. Us- tothebasinofattractionsurroundingaminimumonthe
ing this PES transformation, the non-icosahedral global originalPES.ThisPESisthensearchedbyconstanttem-
minima at N=38, 98,102–104 were much more likely to perature Monte Carlo. Additionally, the algorithm has
be found by their multi-start minimization algorithm.41 been found to be more efficient for clusters if the config-
An additional transformationhad to be applied in order urationisresettothatofthenewlocalminimumateach
to find the global minima at N=75–77. accepted step.46
Itisthereasonsforthesuccessofthissecondapproach There are two ways that one might incorporate a fur-
thatweexamineinthispaper. Inparticularweshowhow ther PES transformation into this algorithm. One could
Schoen and Locatelli’s transformation affects the topog- use basin-hopping to first find the global minimum of
raphy of the PES.We also show how the transformation the transformed PES, then reoptimize the n lowest
low
can be incorporated as an element of an existing algo- energy minima under the original potential. However, if
rithm, namely basin-hopping. theglobalminimumoftheoriginalPESisnotamongthe
n lowest energy minima of the transformed PES this
low
approachis bound to fail.
II. METHODS Alternatively, at each step one could first optimize a
new configuration using the transformed potential, then
reoptimize the resulting minimum using the original po-
The atoms in the clusters interact via the Lennard-
Jones potential:42 tential. By incorporating this second minimization the
shortcomings of the first approachare avoided. Further-
12 6 more, if the energy of this final minimum is used in the
σ σ
E =4ǫ − , (1) Metropolis acceptance criterion, the Boltzmann weight
LJ
r r
i<j"(cid:18) ij(cid:19) (cid:18) ij(cid:19) # ofeachminimumisunchanged. However,theoccupation
X
probabilityofaparticularminimumwillbe proportional
whereǫisthepairwelldepthand21/6σistheequilibrium theareaofthebasinofattractionoftheminimumonthe
pair separation. To this, Schoen and Locatelli added a transformed rather than the original PES, i.e.
2
p ∝n A˜ exp(−βE ), (4) we obtain a good representation of the lower energy re-
i i i i
gions of the PES. At each value of µ we obtained
comp
where ni is the number of permutational isomers of i a sample of 6000 minima. The effect of µcomp on the
and A˜i is the total area of the basins of attraction of number of stationary points, which we noted for CLJ13,
the minima on E˜ which when reoptimized on E lead to is again evident (Table II). As µ increases, n ,
comp search
minimum i. Therefore, if the relative area of the global the number of minima from which we have to perform
minimumislargeronthetransformedPES,optimization transition state searches in order to generate the 6000
shouldbeeasierusingthisapproach. Werefertothisver- minima, increasesandit becomes morelikely thata new
sion of the basin-hopping algorithm as two-phase basin transitionstatedoesnotconnecttoanewminimum,but
hopping. This variation is not much more computation- rather to one already in our sample.
allydemandingthanstandardbasin-hoppingbecausethe
starting point for the second minimization is likely to be
close to a minimum of the untransformed PES. 2
There is one further difference from previous imple-
mentations of the basin-hopping algorithm. Previously,
1
we had performed the minimization in Equation (3) by
conjugate gradient.47 However, we have since found a 31/ 34
limitedmemoryBFGSalgorithmthatismoreefficient.48 ) / N 0 38 147
III. RESULTS -Q compave-1 55 75
Q
(
Inglobaloptimizationthe aimoftransformingthe po- 98
-2 102 135
tentialenergysurfaceistomaketheglobalminimumeas-
iertolocate. Typically,onethereforewantsthetransfor- 13
mationtoreducethe numberofminima andthe barriers -3
20 40 60 80 100 120 140
between them. Furthermore, if the transformation is to N
changetherelativeenergiesoftheminima,onewantsthe FIG. 2. Qcomp for the LJN global minima. To make
the size-dependence more clear the zero is taken to be the
energetic bias towards the global minimum to increase.
As the number of minima and transition states on function, Qave, a four paramater fit to the Qcomp values.
Qave = 25.915N −166.956N2/3 +382.765N1/3 −293.972.
the CLJ PES is small enough that virtually all can
13 Alsoincludedinthefigureareisolated datapoints(crosses)
be found, we canexamine whether the compressiveterm
correspondingtothenon-globalminimaillustratedinFigure
hasthefirstoftheaboveeffectsbyexaminingCLJ asa
13 1 and the second lowest energy minima for N=76, 77 and
functionofµ . Thenumberofminimaandtransition
comp 102–104.
statesclearlydecreasesasµ increases(TableI). Itis
comp
interestingtonotethatminimawithlowsymmetrypref- TheseconddesiredeffectofaPEStransformationisto
erentiallydisappear. The PEStransformationplacesthe change the energetics in a manner that makes the global
cluster in a harmonic potential about its centre of mass. minimummorefavourable. We cangetasimple guideas
Thispotentialplaysarolesimilartoasoftsphericalbox, tohowtheenergiesofthe minima dependonµ ifwe
comp
and so less compact minima disappear from the PES as assume there is no structural relaxation in response to
µ increases. Similar results are found when periodic changing µ . Then E =E +µ Q where
comp comp CLJ LJ comp comp
boundaryconditionsareapplied—thenumberofminima the order parameter, Q = |r − r |2/σ2, is
comp i i c.o.m.
is much less than for a LJ cluster of equivalent size and evaluated at µ =0. From the values of Q we can
comp comp
P
the number of minima decreases as the pressure in the predict the changes in the relative energies of any two
cell is increased.49,50 minima.
Itisalsoworthnotingthatthemagnitudeofthedown- Q is a measure of the compactness of the cluster,
comp
hill barriersrelativeto the energydifference betweenthe andfromFigure2onecanseehowthecompactnessofthe
minima decreases as µ increases (Table I). In the global minima depends on size. For the first two shells
comp
terminology used by Berry and coworkers,51 the pro- the icosahedral global minima are most compact when
files of the pathways to the global minimum become complete Mackay icosahedra can be formed, e.g. N=13
morestaircase-likeandlesssawtooth-likewithincreasing and 55. However, for the third shell the most compact
µ . The combination of the changes to the number icosahedral structure is at N=135, where twelve vertex
comp
of stationary points and the barrier heights act to make atomsoftheMackayicosahedronaremissing,ratherthan
relaxation to the icosahedral global minimum easier as at N=147.
the PES is further transformed. If we examine LJ as an example of a cluster with
38
Next, we examine the CLJ cluster. For a cluster of a non-icosahedral global minimum, we see that this size
38
thissizeitisnotfeasibletoobtainacompleterepresenta- corresponds to a pronounced minimum in Figure 2—the
tion of the PES in terms of stationary points, so instead truncated octahedron is particularly compact compared
3
to the other global minima of similar size. Further- the global minimum to which the other minima directly
more, from Figure 3b we can see that the LJ global join. Foramultiple-funnelPEStherewouldbeanumber
38
minimum has the lowest value of Q of all the LJ of major stems which only join at high energy.
comp 38
minima. Therefore, the energy gap between the global
minimum and the lowest-energy icosahedral minimum From the disconnectivity graph of LJ38 one can de-
increases with µ (Table II). To visualize how this duce that the cluster has a double-funnel PES (Figure
comp
deepening of the fcc funnel changes the PEStopography 4a). There is a narrow funnel associated with the global
we present disconnectivity graphs of CLJ for a range minimum, and a wider funnel associated with the icosa-
38
of µ values in Figure 4. hedral minima. There are a number of low-energy min-
comp
Disconnectivity graphsprovidea representationofthe ima at the bottom of the icosahedral funnel, which, al-
barriers between minima on a PES.52,53 In a disconnec- though they have only small differences in the way the
tivity graph, eachline ends at the energy of a minimum. outerlayerisarranged(e.g.thesecondlowesticosahedral
At a series of equally-spaced energy levels we compute minimum,38C,isdepictedinFigure1),canbeseparated
which(setsof)minimaareconnectedbypathsthatnever by moderate-sizedbarriers. As a resultthere is a certain
exceed that energy. We then join up the lines in the dis- amountoffinestructureatthebottomoftheicosahedral
connectivity graph at the energy level where the corre- funnel with not allminima joined directly to the stem of
sponding (sets of) minima first become connected. In a the lowest-energy icosahedral minimum. From the data
disconnectivity graph an ideal single-funnel PES would inTableIIonecanseethattherearemanymoreminima
berepresentedbyasingledominantstemassociatedwith associated with the icosahedral funnel.
(a) (b)
85 102
84
100
83
82
98
Q 81 Q
comp comp
80 96
79 C
94
78 A
77 H 92
A
-150 -149 -148 -147 -146 -145 -144 -174 -173 -172 -171 -170 -169 -168 -167 -166
Energy /ε Energy /ε
(c) (d)
470
310
468
466
305
464
Qcomp Q 462
300 comp
460
B
458
295
456
A
C
454 A
-398 -396 -394 -392 -390 -388 -386 -384 -382 -544 -542 -540 -538 -536 -534 -532
Energy /ε Energy /ε
FIG.3. ScatterplotsofQcomp againstminimumenergyforlargesamplesofminimafor(a)LJ34,(a)LJ38,(a)LJ75 and
(d)LJ98. TheminimadepictedinFigure1arelabelledbythelettercorrespondingtotheirenergeticrank. In(d)thereare
twosubsets: diamondscorrespondtominimafoundwhenthesearchwasstartedfromthetetrahedralglobalminimumand
crossescorrespondtothesetwhenstartedfromthelowest energyicosahedralminimum. Thereisnooverlapbetweenthese
twosetsbecausenopathwaysconnectingthetwofunnelswere located. Thepatternsofpointsfor Q = r /σ are
linear i<j ij
virtually identical to those of this figure, showing that the current transformation is effectively equivalent to Locatelli and
Schoen’s. P
4
(a) (b)
-167.5 -144.0
-168.0 -144.5
-145.0
-168.5
-145.5
-169.0
-146.0
-169.5
-146.5
-170.0
-147.0
-170.5 -147.5
-148.0
-171.0
-148.5
-171.5
-149.0
-172.0
-149.5
-172.5 icosahedral
-150.0
-173.0 -150.5 fcc
C B
fcc
-173.5 icosahedral -151.0
A -151.5
-174.0
(c) (d)
280
-75
278
-76
276
-77
274
-78
272
-79
icos
270
-80
icosahedral
268
-81
266
-82
fcc
fcc
264
-83
262
-84
FIG.4. Disconnectivity graphs of CLJ38 for µcomp = (a) 0 (b) 0.25ǫ (c) 1.0ǫ and (d) 5ǫ. In (a) the 150 lowest-energy
minima are represented in the graph, and in (b)–(d) the 250 lowest-energy minima are represented. The icosahedral and
fcc funnels are labelled. The unitsof energy on thevertical axis are ǫ.
5
280 global minimum.
1
Given the above, it is unsurprising that two-phase
260
2 basin-hoppingfinds the globalminimum morerapidly as
240 µ increases (Figure 6b). At large µ the first-
0.5 comp comp
220 5 passage time is 40 times shorter than for LJ38. Con-
versely, the first-passage time to reach the icosahedral
200
minimum 38C increases. These changes are driven by
C / k
v180 0.25 changes to A˜i in Equation (4). The basin of attraction
of the global minimum increases in size relative to those
160
of the icosahedral minima as the PES is further trans-
140 formed.
0
120
(a)
100
0 0.1 0.2 0.3 0.4 0.5 0.6
temperature /εk-1 1000 34A
FIG.5. Heat capacity curves for CLJ38 with the values
ofµcomp/ǫ,aslabelled. Thecurveswerecalculatedfromour eps
smaemthpoleds.5o4f,556000 minima using the harmonic superposition me / st
e ti
g
Asthefccfunnelbecomesdeeperwithincreasingµcomp assa 100
itincreasesinsizerelativetotheicosahedralfunnel(Fig- p
ure 4). By µcomp=5ǫ the fcc funnel dominates the PES, first
and the disconnectivity graph has the form expected for
34H
an ideal single funnel with only a very small sub-funnel
for the icosahedral minima. These changes are also re-
10
flected in the number of minima associated with both
0 2 4 µ / ε6 8 10
funnels (Table II). comp
(b)
These changes to the PES topography of course affect 10000
the thermodynamics. For LJ there are two peaks in
38
the heat capacity curve (Figure 5). The first is due to
a transition from the fcc global minimum to the icosa- 38C
hedral minima, which is driven by the greater entropy eps
of the latter. The second corresponds to melting. The me / st1000
first transition hinders global optimization because it is e ti
thermodynamically favourable for the cluster to enter g
a
the icosahedral funnel on cooling from the molten state, ass
p
where it canthen be trapped.37,38 However,as µcomp in- first
creases, the decreasing entropy of the icosahedral funnel
can no longer overcome the increasing energy difference 100 38A
between the global minimum and the icosahedral funnel
(Table II)andsothis firsttransitionis suppressed. Con-
0 2 4 µ / ε6 8 10
sequently,theheatcapacitycurvesfortheCLJ38 clusters comp
inFigure5showonlyonepeak,indicatingthattheglobal FIG.6. The µcomp-dependence of the first-passage time
minimum is most stable up to melting. (in MC steps) to find the specified minima of (a) LJ34 and
(a) LJ38 from a random starting configuration in two-phase
Of course, the changes to the PES topography and basin-hopping runs. Each point represents an average over
thermodynamics mean that onrelaxationdownthe PES 400 runs. The temperature used is 1.0ǫk−1.
the systemis morelikelyto enter the fcc funnel asµ
comp
increases. Furthermore, the energy barrier to escape Locatelli and Schoen’s transformation works for LJ
38
from the icosahedral funnel relative to the energy dif- becausetheglobalminimumisthemostcompactspheri-
ference between the bottoms of the two funnels becomes calminimum. However,this doesnotnecessarilyhaveto
smaller (Table II), thus making escape from the icosahe- be the case, even for those clusters with non-icosahedral
dralfunneleasier. Toquantifytheseeffectsweperformed global minima. From Figure 2 one can see that the non-
annealing56 simulations for CLJ at a number of values icosahedral global minima at N=98 and 102–104 have
38
of µ (Table III). For LJ 80% of the longer an- particularly low values of Q and Figure 3d confirms
comp 38 comp
nealing runs ended atthe bottom of the icosahedralfun- that the Leary tetrahedron, 98A, has the lowest Q
comp
nel, and only 2% at the global minimum. However, by value of all the LJ minima. Therefore, Locatelli and
98
µ =5ǫ 99.5% of the long annealing runs reached the Schoen were able to locate these global minima. How-
comp
6
ever, for N=75-77 the values of Q for the Marks global minimum is only one of the more compact min-
comp
decahedra are not set apart from the nearby icosahedral ima, the transformation is less beneficial for global opti-
globalminima (Figure2)andFigure3cshowsthatthere mization. By contrast, for sizes with icosahedral global
are a number of LJ minima which have lower values of minima the transition is often unhelpful, as we saw for
75
Q than75A.Inparticular,the icosahedralminimum LJ , because the global minimum is much less likely to
comp 34
75C that is third lowest in energy has a lower Q , be the most compact structure. Therefore, the trans-
comp
and the Marks decahedronis no longerthe CLJ global formation needs to be used in combination with other
75
minimum beyond µ =3.1ǫ. methods. As the transformationis mostlikely to be suc-
comp
ThegeometricrootofthisbehaviouristhattheMarks cessful for clusters where other methods fail it can act
decahedraatN=75–77aretheleastsphericalofthenon- as a good complement to them. For example, when the
icosahedral global minima. The 75-atom Marks deca- basin-hopping algorithm is applied, usually a series of
hedron is somewhat oblate and some of the icosahedral runs are performed at each size. If one of the runs used
minima with which 75A is competing are prolate by a the two-phase approach, this would increase the chance
similar degree, leading to comparable values of Q . of success for those sizes where the PES had a multiple-
comp
Therefore, although the transformation may aid global funnel topography.
optimization by reducing the number of minima and by Other PES transformations could also be usefully em-
increasing the energy of many minima relative to the ployedalongsidestandardbasin-hoppingrunsinthistwo-
Marksdecahedron,unlikeforLJ itdoesnotremovethe phase approach,if they are likely to aidglobaloptimiza-
38
fundamental double-funnel character of the PES. To lo- tion for some sizes. For example, increasing the range
catetheglobalminimumLocatelliandSchoenhadtoadd of the potential is another transformation that reduces
an additional ‘diameter penalization’ to the potential.41 the numberofstationarypointsonthe PES.43 Usingthe
Locatelli and Schoen found that for many of the clus- transformations alongside standard runs avoids one of
terstheirtransformationdidnotaidglobaloptimization. the major difficulties associated with PES transforma-
Thiswasnotunexpected,butsimplyreflectsthefactthat tions. They are rarely universally effective, but rather
often the icosahedral global minima are not the most there are likely to be some instances when they destabi-
compact minima. We analyse one example. At N=34 it lizetheglobalminimum,thusmakingoptimizationmore
is possible to forma compact Learytetrahedron(34H in difficult. This is certainly the case when increasing the
Figure 1), which is the eighth lowest-energy LJ mini- rangeofthepotential,wheretherange-dependenceofthe
34
mum. This structure has a significantly lower value of most stable cluster structure is well-documented.10,57
Q than the global minimum (Figure 3a). As a re- Although we have seen how a compressive transfor-
comp
sult, the Leary tetrahedron becomes the CLJ global mation can be useful in aiding the global optimization
34
minimumatµ =0.3ǫ. Theresultsoftwo-phasebasin- of LJ clusters, an important question is how generally
comp
hoppingrunsaresimilartothoseforLJ inthatasµ useful it will be. Although this question can only be
38 comp
increasesthecompactnon-icosahedralstructurebecomes definitively answered through applications to a variety
significantly easier to locate and the low-energy icosahe- of systems, one would expect it to be useful for metal
dral minima more difficult (Figure 6). The difference, and simple molecular clusters that form compact struc-
though, is that now this scenario is undesirable, because tures,particularlythosethatfavour12-coordination. For
it is the global minimum that is becoming more difficult these systems, as with LJ clusters, the strength of this
to reach. approachwouldbelocatingthoseglobalminimathatare
not based on the dominant morphology, because the al-
ternative morphologies are only likely to be most stable
IV. CONCLUSIONS when they are compact and sherical. It might also be
usefulinsystemssuchasproteinswheretherearealarge
number of less compact unfolded configurations. How-
By analysing the effect of a compressive transforma- ever, it would not be useful for clusters of substances,
tion on the PES topography we have obtained insights suchaswaterandsilicon,whichformopennetworkstruc-
intothereasonsforitssuccessinaidingthe optimization tures where the liquid can be denser than the solid.
of LJ clusters that have non-icosahedral global minima.
Firstly, we have shown that the transformation reduces
the number of minima and transitionstates on the PES.
ACKNOWLEDGMENTS
Secondly, for examples where, as is often the case, the
non-icosahedral global minimum is the most compact
structure, the transformation causes the funnel of the J.P.K.D. is the Sir Alan Wilson Research Fellow at
global minimum to become increasingly dominant. For Emmanuel College, Cambridge. The author is grateful
LJ the PES has a double funnel, whilst at large µ to David Wales for supplying a modified version of the
38 comp
the PES has an ideal single-funnel topography, enabling basin-hopping code, and would also like to thank Marco
the system to relax easily down the PES to the global LocatelliandFabioSchoenforhelpfuldiscussionsandfor
minima. However, when as for LJ , the decahedral sharing results prior to publication.
75
7
31S.K.Gregurick,M.H.Alexander,andB.Hartke,J.Chem.
Phys.104, 2684 (1996).
32J. A. Niesse and H. R. Mayne, J. Chem. Phys. 105, 4700
(1996).
33W. Pullan, Comp. Phys.Comm. 107, 137 (1997).
1D.J.WalesandH.A.Scheraga,Science285,1368(1999). 34M.D.WolfandU.Landman,J.Phys.Chem.A102,6129
2F. H. Stillinger, Phys.Rev.E 59, 48 (1999). (1998).
3L. T. Wille, in Annual Reviews of Computational Physics 35C. Barr´on, S. G´omez, D. Romero, and A. Saavedra, Appl.
VII, edited by D. Stauffer (World Scientific, Singapore, Math. Lett. 12, 85 (1999).
2000). 36B. Hartke,J. Comp. Chem. 20, 1752 (1999).
4L. T. Wille, Chem. Phys. Lett. 133, 405 (1987). 37J. P. K. Doye and D. J. Wales, Phys. Rev. Lett. 80, 1357
5J. A. Northby,J. Chem. Phys. 87, 6166 (1987). (1998).
6T. Coleman and D. Shalloway, J. Global Optimization 4, 38J. P. K. Doye, D. J. Wales, and M. A. Miller, J. Chem.
171 (1994). Phys.109, 8143 (1998).
7G. L. Xue,J.Global Optimization 4, 425 (1994). 39D.J. Wales, personal communication.
8S. Gomez and D. Romero, Proceedings of the First Euro- 40B. Hartke,personal communication.
pean Congress of Mathematics (Birkhauser, Basel, 1994), 41M. Locatelli and F. Schoen, Computational Optimization
Vol. III,pp.503–509. and Applications submitted (1999).
9J.Pillardy andL.Piela, J.Phys.Chem. 99,11805 (1995). 42J. E. Jones and A. E. Ingham, Proc. R. Soc. A 107, 636
10J.P.K.Doye,D.J.Wales,andR.S.Berry,J.Chem.Phys. (1925).
103, 4234 (1995). 43M. A. Miller, J. P. K. Doye, and D. J. Wales, J. Chem.
11J.P.K.DoyeandD.J.Wales,Chem.Phys.Lett.247,339 Phys.110, 328 (1999).
(1995). 44C. J. Cerjan and W. H. Miller, J. Chem. Phys. 75, 2800
12D.M.Deaven,N.Tit, J. R.Morris, and K.M.Ho, Chem. (1981).
Phys.Lett. 256, 195 (1996). 45Apublicdomainversionofthebasin-hoppingcodeisavail-
13D.J.WalesandJ.P.K.Doye,J.Phys.Chem.A101,5111 able from http://brian.ch.cam.ac.uk/software.html.
(1997). 46R.P.WhiteandH.R.Mayne,Chem.Phys.Lett.289,463
14C. Barr´on, S. G´omez, and D. Romero, Appl. Math. Lett. (1998).
10, 25 (1997). 47W. H. Press, B. P. Flannery, S. A. Teukolsky, and
15R.H. Leary, J. Global Optimization 11, 35 (1997). W. T. Vetterling, Numerical Recipes (Cambridge Univer-
16D.Romero,C.Barr´on,andS.G´omez,Comp.Phys.Comm. sity Press, Cambridge, 1986).
in press (1999). 48D.Liu and J. Nocedal, Mathematical Programming B45,
17R. H. Leary and J. P. K. Doye, Phys. Rev. E 60, R6320 503 (1989).
(1999). 49T.A.WeberandF.H.Stillinger,J.Chem.Phys.80,2742
18D. J. Wales, J. P. K. Doye, A. Dullweber and F. (1984).
Y. Naumkin, The Cambridge Cluster Database, URL 50A.Heuer, Phys.Rev.Lett. 78, 4051 (1997).
http://brian.ch.cam.ac.uk/CCD.html. 51K.D. Ball et al.,Science 271, 963 (1996).
19URL 52O. M. Becker and M. Karplus, J. Chem. Phys. 106, 1495
http://www.vetl.uh.edu/∼cbarron/LJ cluster/researchpot.html. (1997).
20M.R.HoareandJ.McInnes,FaradayDiscuss.,Chem.Soc. 53D. J. Wales, M. A. Miller, and T. R. Walsh, Nature 394,
61, 12 (1976). 758 (1998).
21C. J. Tsai and K. D. Jordan, J. Chem. Phys. 99, 6957 54D.J. Wales, Mol. Phys. 78, 151 (1993).
(1993). 55G. Franke, E. R. Hilf, and P. Borrmann, J. Chem. Phys.
22J. P. K. Doye, M. A. Miller, and D. J. Wales, J. Chem. 98, 3496 (1993).
Phys.111, 8417 (1999). 56S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, Science
23P. E. Leopold, M. Montal, and J. N. Onuchic, Proc. Natl. 220, 671 (1983).
Acad. Sci. USA89, 8271 (1992). 57J. P. K. Doye and D. J. Wales, J. Chem. Soc., Faraday
24J. D. Bryngelson, J. N. Onuchic, N. D. Socci, and P. G. Trans. 93, 4233 (1997).
Wolynes, Proteins: Structure, Function and Genetics 21,
167 (1995).
25A.L. Mackay, ActaCryst. 15, 916 (1962).
26L. D.Marks, Phil. Mag. A 49, 81 (1984).
27J. P. K. Doye, M. A. Miller, and D. J. Wales, J. Chem.
Phys.110, 6896 (1999).
28M. A. Miller, J. P. K. Doye, and D. J. Wales, Phys. Rev.
E 60, 3701 (1999).
29Z. Li and H. A. Scheraga, Proc. Natl. Acad. Sci. USA 84,
6611 (1987).
30D. M. Deaven and K. M. Ho, Phys. Rev. Lett. 75, 288
(1995).
8
TABLE I. The number of minima, nmin, and transi- TABLEII. Properties of the CLJ38 PES for a sample of
tion states, nts, for CLJ13 as a function of µcomp. For 6000 connected minima as a function of µcomp. nts is the
each minimum 30 transition state searches were performed; numberoftransitionstatesconnectingtheseminima. ∆E is
these searches were parallel and antiparallel to the eigen- the energy difference between the global minimum and the
vectors with the fifteen lowest eigenvalues. ∆E, bu, bd lowestenergyicosahedralminimumandbfcc (bicos)istheen-
are the average energy difference, uphill barrier and down- ergy barrier that has to be overcome to escape from the fcc
hill barrier, respectively, where the average is over all the (icosahedral) funnel and enter the icosahedral (fcc) funnel.
non-degeneraterearrangementpathways. (Degeneratepath- Of course, ∆E = bfcc−bicos. For the nsearch lowest-energy
ways connect different permutational isomers of the same minima 20 transition state searches were performed; these
minimum.) ∆E =b −b . searches were parallel and antiparallel to the eigenvectors
u d
µcomp/ǫ 0 0.5 1 2.5 5 10 25 wbeitrhs otfhemtinenimloawiensttheeigfecncvaanludesic.onsafchceadnradlnfuicnosnealrseatthtehneuemn--
nmin 1467 769 470 169 75 33 10 ergy at which the two funnelsbecome connected.
nts 12435 5820 3010 801 262 100 37
∆E/ǫ 1.593 3.172 4.501 7.191 11.215 20.701 40.176 µcomp/ǫ 0 0.25 0.5 1 2.5 5
bu/ǫ 2.201 3.939 5.396 8.231 12.346 21.759 42.263 nts 8633 9111 9911 11656 17137 23270
bd/ǫ 0.609 0.767 0.896 1.041 1.131 1.058 2.087 nsearch 1271 1277 1491 1924 3107 4253
b /∆E 0.382 0.242 0.199 0.145 0.101 0.051 0.052 ∆E 0.676 1.550 2.274 3.564 6.120 9.893
d
b /ǫ 4.219 4.795 5.256 6.143 8.892 12.659
fcc
bicos/ǫ 3.543 3.245 2.981 2.580 2.772 2.766
bicos/∆E 9.893 2.094 1.311 0.724 0.453 0.280
n 92 113 73 106 104 86
fcc
nicos 912 439 194 27 5 6
nfcc/nicos 0.11 0.26 0.38 3.93 20.8 14.33
TABLEIII. ResultsofannealingsimulationsforCLJ38as
afunctionofµcomp. fOh(ncycles)isthefractionoftheanneal-
ingrunsthatterminatedattheglobalminimum,andficos is
thefraction ofrunsthatendedinthelowestfiveicosahedral
minimum. Each annealing run involves a linear decrease in
thetemperaturefromtheliquidto0Kinn MonteCarlo
cycles
cycles. The results are averages over 200 annealing runs.
µcomp/ǫ 0 0.25 0.5 1 2.5 5
f (106) 0% 2.5% 7.5% 19.5% 67% 79.5%
Oh
f (107) 2% 14% 31% 66.5% 97% 99.5%
Oh
ficos(106) 37% 29.5% 12.5% 7.5% 1% 0%
ficos(107) 80% 56.5% 38% 6.5% 0% 0%
9