Genetic operators for the atomic cluster problem

Computer Physics Communications ELSEVIER

Computer Physics Communications !07 (1997) ! 37-148

Genetic operators for the atomic cluster problem W.J. P u l l a n Department of Mathematics and Computing, Central Queensland University. Rockhampton. Queensland. Australia. 4702 Received 7 April 1997

Abstract This paper investigates the recent advances in applying genetic algorithms to non-bonded molecular conformation problems. Through the use of deterministic local optimisation to reduce the search domain and innovative phenotype genetic crossover operators, genetic algorithms are now able to efficiently solve problems such as detennining the minimum energy configurations of pure atomic clusters (N = 2 . . . . . 100), mixed atomic clusters (N = 2 . . . . . 20) and molecular clusters (N = 2 . . . . . 32). A comparison of the more conventional genotype crossover operators with the more recently developed phenotype operators is presented. @ 1997 Elsev~,erScience B.V.

PACS: 36.40.-c Keywords: Global optimisation;Genetic algorithm; Molecularconformation

1. Introduction Recently, genetic algorithms have been successfully applied to molecular conformation problems involving clusters of carbon atoms [1], pure atomic clus. ters [2], mixed atomic clusters [3], and clusters of benzene and water molecules [ 4 ]. These genetic algorithms were innovative in that they restricted the search domain through the use of local (deterministic) optimisers and also utilised innovative genetic crossover operators operating in the problem (phenotype) domain rather than the more usual parameter (genotype) domain. The primary objective of this paper is to evaluate these phenotype crossover operators and compare them with the more standard genotype crossover operators, when applied to atomic and molecular cluster problems. Section 2 of this paper overviews the atomic cluster optimisation problem used as a test case while Section 3 describes the genetic algorithm used in this

study. Section 4 presents the results obtained using both genotype and phenotype crossover operators and also an analysis of these two crossover operators. Finally, Section 5 contains a review, future directions and the conclusion.

2. Atomic clusters Atomic clusters are currently an active field of theoretical and experimental research. Atomic clusters are aggregates of atoms, sufficiently small so that a significant proportion of the atoms are present on the surface of the cluster. Theoretical investigations of atomic clusters address the following optimisation problem:

0010-4655/97/$17.00 (~) 1997 Elsevier Science B.V. All rights re~r~.,~'~ Pll SO0 i 0-4655 ( 97 ) 00092- I

Given N particles, interacting with two-body central forces, find the configuration in threedimensional Euclidean space for which the total potential energy attains its global minimum.

W.J. Pullan/Computer Physics Communications 107 (1997) 137-148

138

1

-1

o'.~

0

;

,'5 f

~

2's

Fig. i. Scaled Lennard-Jonespair potential. The potential energy (V) of an N atom cluster can be written as N-i

N

i=i j=i+ l

where r 0 is the Euclidean distance between atoms i and j and v ( r i j ) is the pair potential acting between atoms i and j. For this study, the scaled Lennard-Jones pair potential, v ( r ) = r -12 - 2r -6, having the form shown in Fig. 1, was utilised. Prior to 1987, the most extensive study of atomic clusters was the work of Hoare et al. [ 5-8 ] who developed a general growth algorithm and used it to generate large numbers of stable structures, mainly for N < 55. These were compared to find the lowest energy structures which, in turn, became candidates for the absolute minimal structures. Hoare and Pal [7] observed that, while what they termed as the "icosahedral growth sequence" did not, in general, produce minimal structures, icos~hedral subunits did appear regularly in relaxed configurations generated by other sequences. The icosabedrai lattice can be described as 20 slightly flattened tetrahedrally shaped face-centredcubic units with 12 vertices on a sphere centred at the origin. For the icosahedral lattice, the total number of lattice points on each layer is 1, 12, 42, 92, .... Therefore, the number of lattice points in the sequence of closed shell icosahedral lattices is 1, 13, 55, 147, .... For smaller Lennard-Jones clusters, globally minimum configurations are

• two atoms - separated by one unit; ® three atoms - an equilateral triangle; ® four atoms - a regular tetrahedron; • five atoms - a triagonal bipyramid, slightly contracted along the symmetry axis and distended in the symmetry plane; • six atoms - the regular octahedron with slightly contracted sides; • seven atoms - a pentagonal bipyramid. For more than seven atoms, minimum energy structures are generally based on icosahedral lattice structures as shown in Fig. 2. As an example of the optimal structure of larger atomic clusters, Fig. 3 shows the optimal structure for the 55 atom cluster. In this structure, one atom is positioned at the centre, 12 are positioned on the vertices of an icosahedral lattice which is centred on the first atom and the remaining 42 are positioned at the vertices of a second enclosing icosahedral lattice which is also centred on the first atom. From [ 5], the total number of distinct local minima for N = 2 . . . . . 13 atom clusters, using a scaled Lennard-Jones pair potential, is 1, I, 1, 1, 2, 4, 8, 18, 57, 145, 366, and 988. As stated in [5], this data is approximated by the equation re(N) = exp(-2.5176 + 0.3572N + 0.0286N 2) which, if correct, predicts that the potential energy surfacc for a 100 atom cluster has of the order of 10140 local minima! A number of optimisation methods have been applied to the atomic cluster problem. These include Lattice Optimisation/Relaxa~ion Techniques [ 9,10], Simulated Annealing [11,12], DC Transformation/GOP Algorithm [ 13], and the Diffusion Equation Method [ 14,15 ]. Of these, the Lattice Optimisation/Relaxation techniques have been most successful and currently provide most of the definitive results. Of particular relevance to this study, Hartke [ 16] used a genetic algorithm to find global minimum energy structures for four atoms. An encoding method, using just four parameters to specify the four atom configuration, allowed a standard GA to find the global minimum energy. Mestres and Scuseria [ 17] also used a GA to find global minimum energy configurations for Lennard-Jones atomic clusters of up to 13 atoms and Gregurick et al. [ 18] employed a modified deterministic/stochastic genetic algorithm to find all currently accepted global minima for Lennard-Jones

W.J. Pullan/Computer Physics Communications 107 (1997) 137-148

(a) N = 8, e n e r g y = - 19.822.

139

(b) N = 9, e n e r g y = - 2 4 . 1 1 3 . f

(c) N = 10, e n e r g y = - 2 8 . 4 2 3 .

(d) N = 11, e n e r g y = - 3 2 . 7 6 6 .

(e) N = 12, e n e r g y = - 3 7 . 9 6 8 .

(f) N = 13, e n e r g y = - - 4 4 . 3 2 7 .

Fig. 2. Globally optimal structures for N = 8 . . . . . 13 atom scaled Lennard-Jones clusters.

atomic clusters in the range 2-29. However, their technique did require a seed growth method for generating initial structures for all cluz~ers containing more than 20 atoms. More recently, Deaven et al. [ 2 ], using a real-coded GA with innovative phenotype ge-

netic crossover operators, were able to find all global minimum energies for Lennarddones clusters in the range 2-100 including new minima for 38, 65, 69, 76, 88 and 98 atom clusters.

140

W.J. Pullan/ Computer Physics Communications 107 (1997) 137-148

Fig. 3. Optimal structure for the 55 atom scaled Lennard-Joaes cluster. While this is a non-bonded structure, bonds have been included in the figure to highlight the basic icosahedral lattice structure of the optimal cluster.

3. Genetic Algorithms Genetic Algorithms (GAs) are based on the process of natural selection where the proportion of "fitter" individuals in a population tends to increase at each generation [ 19]. Standard GAs operate on a binary encoded version of the parameters of a problem and use genetic operators such as proportional selection (based on a fitness factor), random mutation and crossover to generate a new population from an existing population of individuals. Theoretical analysis and experimental results [19] show that standard (or binary encoded) GAs are efficient when applied to problems whose parameters can be encoded as short, low-order schema which are relevant to the underlying problem and relatively unrelated to schema over other fixed positions [ 19 ]. ltowever, the atomic cluster optimisation problem meets neither of these criteria because ( I ) Atomic positions are real valued (floating point) parameters. Using a binary coded GA requires that these positions be discretised so they can be coded as binary integers. To obtain an accuracy of 3 digits after the decimal place, for a 30 atom cluster, requires a relatively large binary solution vector of length 1260 bits. This in turn generates a search space of approximately 1042oelements. (2) In evaluating the potential energy (and hence the fitness factor) of a water cluster, the position of each molecule is only meaningful in the context of the positions of all other atoms. That is, schema over different positions are not unre-

lated, so standard GAs can be expected to perform poorly when applied to the molecular cluster problem [ 19]. As the theoretical basis for GAs was originally restricted to standard GAs, the alternative of using GAs whose parameters were encoded as real numbers has been relatively infrequently chosen. However, it has been shown recently [20] that the "implicit parallelism" [ 19 ] of standard GAs does not depend on binary encoding and there is a growing body of evidence that GAs with real parameters have acceptable performance [21,22] provided that suitable genetic operators are utilised. With regard to molecular structure problems, this was first demonstrated by Deaven and Ho [ l ] where they addressed the first problem listed above by the use of a local (deterministic) optimiser which had the effect of mapping a range of parameter values (those corresponding to energies within a "catchment basin" on the potential energy surface) to a particular set of parameter values (that corresponding to the local minimum of the catchment basin). In effect, the !ocal optimiser performs a natural discretisation cf the parameter values where the interval of discretisation is directly related to the size of a catchment basin. Deaven and Ho addressed the second problem identified above by the use of innovative geometric genetic crossover operators operating in the problem (phenotype) space rather than the more normal genetic crossover operator which operates in the parameter (genotype) space. For molecular structure problems in particular, these geometric crossover operators tend to preserve the structural relationships that have developed within individuals and provide the building blocks necessary to create fitter individuals. Direct evidence of the effectiveness of this type of GA, as compared to the standard GA, is provided by the results obtained for pure atomic clusters in [ 18] as compared to those obtained in [2]. In [! 8], using randomly generated initial clusters with a standard GA, all clusters in the range N = 2 . . . . . 20 were successfully optimised. In contrast, Deaven and Ho [2], using randomly generated initial clusters, were able to obtain all currently accepted global minima for clusters in the range N = 2 , . . . , i00. The effectiveness of the local optimiser in reducing the search space is demonstrated in Fig. 4 where the natural logarithms of 236N (the search space size for a

W.J. Pullan/ Computer Physics Communications 107 (1997) 137-148

141

These three variants were used with equal probability during the optimisation process. The hemisphere variant performed a relatively global search while the quadrant and orthant variants searched closer to the structural vicinity of the parent clusters.

2500

2000

8 ~'1500

3.3. Genetic mutation .~'1000

m

I

Z

500

!

10

20

30

40

50

60

70

80

90

100

Number of atoms sn cluster

Fig. 4. Natural logarithmsof 2?'6N(upper curve), the search .~pace size for a standard binary, encoded GA where each atomic coordinate requires 12 bits and exp( -2.5176 + 0.3572N + 0.0286N2) (lower curve), the search space size when only the local minima are considered. standard binary encoded GA where each atomic coordinate requires 12 bits) and exp( - 2 . 5 1 7 6 + 0 . 3 5 7 2 N + 0.0286N 2) [ 5] (the search space size when only the local minima are considered) are plotted. The GA used to obtain the results presented in this paper contains many of the features of the Deavon and Ho GA [ 1,21 and is now described. 3. !. Genetic encoding All atoms were encoded using three real parameters to specify the Cartesian coordinates of each the atom.

The following equally likely mutation operations were performed, with probability Pm: ® a combination of BFGS local minimisations and APSE [ 23] searches were used to perform a multistart global optimisation within the structural vicinity of a cluster; o a randomly selected atom was removed, a new atom added at a randomly selected position and a subsequent BFGS local minimisation was performed; ® the atom which had the highest binding energy was removed, a new atom added at a randomly selected position and a subsequent BFGS local minimisation was performed. All child clusters that were not mutated were locally optimised using a BFGS optimiser.

3.4. Genetic parameters The population size was always set at twelve with a probability of mutation (Pro) of 0. i and all possible crossovers were performed. No individual was added to the genetic pool if its energy was within 0.0005 of an individual already in the genetic pool.

3.2. Genetic crossover

4. Results The genetic crossover operator for generating new clusters translated each parent cluster so that the centre of mass of the cluster was positioned at the origin of the coordinate system and then randomly rotated the clusters around each axis of the coordinate system. Child clusters were constructed by • using opposing hemispheres (above and below the xy-plane) from each parent cluster; ® using a quadrant from one parent cluster and the remainder of the sphere from the other parent cluster; ® using an orthant from one parent cluster and the remainder of the sphere from the other parent cluster.

In the following subsections, the effectiveness of the standard single-point genotype and the phenotype crossover operators is compared using experimental data. A measure of the effectiveness of crossover operators is then proposed and analytical expressions derived for this measure for both types of crossover operators. These analytical expressions are then verified from experimental data. Finally, based on these observations, a modification to the phenotype crossover is proposed and experimental results showing the resultant improvement presented.

W.J. Pullan/Co,nputer Physics Communications 107 (1997) 137-148

142

r." -20 ® -20

:.

!"

iIi;•

• ~

;."

E -40

i -oL- .

i.

!

." :

"

,'-

,: i . o•

%

$

$

.'.. " i

"i'

?-'" "

".-:

.:•

°

"

.

"

1 -aot-

,"

_~-1~

-12(

_

0~,.

i

-140

- 120

~

I - 100

I

QI

-80 .-600 Average parent cluster energy

t -40

/

I -20

-14( -140

'

:

,

,

,

i

-100

I

-80

i

-60

-

40

i

-20

Average parent cluster energy

(b) X1 - energy after B F G S optimisation.

(a) X z - energy after crossover. ;

i

-120

,

.

-2O

s

...o

-2O

:i -40

-40

(n -60

-60

-80

~ 40

f, ~ -i0c, o

-120

/

"14f~0

'i20

'100

-80

-60

-40

-2()

0

Average parent cluster energy

(c) X2 - energy after crossover.

-140 / -140

J -1:~0

-100

-8~0

-610

-4~0

-210

Average parent cluster energy

(d) X2 - energy after B F G S optimisation.

Fig. 5. Average patent, crossover and BFGS optimised cluster energies for the Xl and X2 crossover operators during the optimisation of an N = 30 atomic cluster. All positive energy clusters have been omitted from this plot. Each point represents a single crossover operation or BFGS oplimisation and shows the relationship between the average parent cluster energy and either the child energy after crossover or the final optimised energy of the child cluster. Points below the solid diagonal line show where the child cluster had a lower energy than the average parent cluster energy.

4.1. Experimental observations To determine the effect of the phenotype genetic crossover operation, the GA was modified so that at each crossover operation, two child clusters were first created using a standard single-point genotype crossover operator (Xl) and then regenerated using the phenotype crossover operator (X2) described in Section 3.2. At the completion of each crossover the

following data was logged for both the X1 and X2 operators (during the optimisation of an N = 30 atomic cluster): ® the average energy of the two parent clusters; , the energies of both child clusters after crossover; ® the BFGS optimised energy of each child cluster. Fig. 5 shows both the distribution of cluster energy immediately after crossover and the BFGS optimised energy of child clusters, as a function of the average

w.z Pullan/Computer Physics Communications 107 (1997) 137-148

energy of the two parent clusters, for the g~ and X2 crossover operators. In all cases, data relating to clusters with positive energy has been omitted. Each point represents a single crossover operation or BFGS optimisation and shows the relationship between the average parent cluster energy and either the child energy after crossover or the final BFGS optimised energy of the child cluster. From the rdative distribution of these points it is clear that the ,1'2 crossover operator is more effective at exploring the search domain at lower energies than the XI crossover. In fact, for X~, no child cluster had an energy lower than the average parent cluster energy immediately after crossover while 13% had an energy lower than the average parent cluster energy after the subsequent BFGS optimisation. For X~, however, 5% of child clusters had an energy lower than the average parent cluster energy immediately after crossover while 21% had an energy lower than the average parent cluster energy after the subsequent BFGS optimisation.

4.2. Crossover analysis

For the atomic cluster problem, one measure of the effectiveness of a genetic crossover operator is the energy of the child clusters as compared to the energy of the two parent clusters. A crossover operator that tended to generate child clusters with high energy (relative to the parent cluster energies) could be viewed as less efficient than one that generated lower energy clusters as it is probable that more computational effort is required in the subsequent BFGS local optimisation to reduce this higher energy. For the purposes of this analysis, however, a simple count of the pairs of atoms which are "close" (within 0.5) in child clusters is used as a measure of the effectiveness of a genetic crossover operator. As shown in Fig. 6a, there is a close relationship between the average cluster energy of clusters and the number of close atoms in the cluster. Clearly higher energy values are, on average, associated with highel" counts of close atoms. As evidence of this measure as an arbiter of success, Fig. 6b shows the relationship between the number of close atoms in a cluster and the average final energy of the cluster after applying the BFGS local optimiser. Clearly an initially low number of close atoms tends to result in a lower final optimised energy.

143

Expressions for the expected number of close atoms after gl and X2 crossovers are now analytically derived using the assumptions that all parent clusters contain N uniformly distributed atoms within an enclosing sphere of radius R with volume V = 47rR3/3 and that no atom in a parent cluster is close to any other atom. In this context, atoms are defined as being close if the distance between the centres of the two ato~as is less than some value r (r -- 0.5 in this study), b'or the purposes of analysis of the crossover operator, it is convenient to view this in terms of the atoms already in the cluster having a "sphere of influence" of radius r while the atom being added is represented by a single point. Two atoms are defined as being close if the point representing the atom being added lies within a sphere of influence of an existing atom. 4.2.1. XI crossover

Given two parent clusters X! and X2, let a child cluster X be created by initially randomly selecting K atoms from parent XI and then incrementally adding N - K atoms from parent X2. The probability Pi that the first atom added from )(2 lies within the sphere of influence of one of the K existing atoms in X is approximately the ratio of the volume occupied by the spheres of influence of the K atoms to the volume into which the new atom may be added, that is, P| = Kr 3/ R 3 .

As no atom within a parent cluster may be close to another atom, the volume into which the second new atom may be added is 41r(R 3 - r 3)/3. However, as the first atom added may have been close to one of the existing K atoms, to a first approximation the expected value of the volume which will cause an overlap is 4~'( Kr 3 - Pt (r 3/2) )/3. Therefore, the probability P2 that the second atom added is within the sphere of influence of one of the existing K atoms is P2 = ( K r3 - Pl( r 3 / 2 ) ) / ( R3 - r3) -

For the third atom, the volume into which it may be added is 4!7"(R 3 - 2r 3 )/3 and the expected value of the volume which will cause an overlap is 4~r(Kr 3 (r 3/2) (,°2+P!) )/3. Therefore, the probability P3 that the third atom added is within the sphere of m,uence "-'~ of one of the existing K atoms is approximately P3 = (Kr3 - (r3/2) (Pi + P2) ) / ( R3 - 2r3) •

WJ. Pulhm/Computer Physics Communications 107 (1997) 137-148

144 X |04

; m II I

,I m

I

ltt

-40

,

•

z

"

•

II II

g s

•

.

3

"

:

;

,

|

•

•

I!

~

-80

8

I

II

I

•

i

,

,

,

!

=

.

.

:

,

,

.

,

•

•

•

:

•

!. / i . .. !; ~: 'i

.

i

•

!

;

.

~

"

: .

,<

"~I

:

"6

ltl

2

.

~

" i, .

-6o

III

g

.

. i

N

;

.

-20

II

Ill Ii

II It

I

,.

I

"

.

I

!i.

-100

|

!

1111111

-120

o

10

5

1

-

_1401

I 2

, 4

Numbet" of close atoms in cluster

i 6

t 8

, 10

t 12

I 14

1

1

20

Number of c;ose atoms in mmal clusler

(a) Average cluster energy.

(b) Distribution of B F G S optimised energy.

Fig. 6. (a) The relationship between the average cluster energy and the number of close atoms (within 0.5) sampled during the complete optimisation of an N = 30 atomic cluster. The correlation coefficient for this data is 0.9476 giving a positive correlation between the average cluster energy and the number of close atoms• (b) The effect that the number of close atoms in an N = 30 atomic cluster has oa the final energy found by ~he BFGS local optimiser. Generally, the higher the number of close atoms, the lower the probability of the local optimiser finding a relatively lower energy. The correlation coefficient between the number of close ,~toms and the average BFGS optimised energy (flagged by * in graph) is 0.9373.

Continuing this process, the probability that atom i is within the sphere of influence of one of the existing K atoms is i-I

7i

......

, . . . . . . . . . . . . T. . . . . . . . .

• .........

.. . . . . . . . . . .

,. . . . . .

~ .........

• .................,

6

gsri

P~ : (K"3 - (r3/2) E P J ) / ( R 3

-

(i-

I)r3).

j=l

Therefore, after N - K atoms have been added, the expected number of close atoms in X is

--,

N-.K

i=1

N-K

= E

i=l

2

gr3

--

i-I ) j=l

(r3/2) ~ P J

3

4

5

6 R

7

8

9

I0

Fig. 7. The expected number of close atoms after a ,l'! crossover as a function of the radius of the minimal sphere which encloses all atoms in the cluster ( N = 40, K = 20).

....... (R 3 - ( i - l)r 3) 4.2.2•

This expression for E[ Nr ] is plotted in Fig. 7 tbr N = 40, K = 20, and varying values of R. It can be seen that the XI crossover operator generates a relatively high number of close atoms, particularly as R decreases (the optimal configuration for the 40 atom cluster can be enclosed by a sphere of radius slightly greater than 2 units).

X2 crossover

Given two parent clusters X! and X2, we need to determine the expected number of close atoms in a child cluster when a volume element Ve is removed from Xi and replaced by the corresponding volume element from X2. For the X2 crossover operator, these volume elements are restricted to hemispheres, quadrants and octants.

W.J. Pullan/Computer Physics Communications 107 ~¢ I.Q 07) !.... ~7- t a8

Given the conditions imposed above on parent clusters, only those atoms in X! and X2 which are within r of a planar face of ~ will be able to be close to another atom in the child cluster X. For an atom within ~ , distance d < r from a planar face and not within r of any other planar face of Ve, the volume of its sphere of influence which protrudes through this planar face is given by

~45

0.5

045t 0.4~

r

Vi, = /

"rr(r 2 -- x 2 ) d x

oos~ ~

o '

d

2

= 27rr3/3 - 1rd(r 2 - d 2 / 3 ) .

wherel <_ i <_ Nt, .

Ignoring corner effects, the expected value for the total volume of all spheres of influence of atoms which protrude through this planar lace of V~ is - di(r 2 - d2/3)).

i=l

If the volume of a cylinder of height r and based on the planar face is denoted by Vc, then the expected number of close atoms in the child cluster is E[ N,.] = N ( V ~ . / V ) ( V s / V c )

N /v

= ( 3 N / 4 R 3) Z ( 2 r 3 / 3

- di(r 2 - d2/3)).

i=!

The expected value c:f N t, is clearly dependent on the shape of Ve and, in particular, the area of the planar contact face. For a hemisphere volume element, Ve = 2zrR3/3 and there is only one planar face, of area 7rR2, .... The ~,,,,~,.,,.a that needs to be consld • ...,,,a -..-v. . . . . . . ., .., ,., ,.~., ~. ,.. ,.. . , - , f atoms within r of this planar face is r

E[ N r l = ( N / V )

B

-7 ~-- ;

Fig. 8. The expected number of close atoms after a X2 crossover as a function of the radius of the minimal sphere which encloses all atoms in the cluster (N = 40, K = 20).

= N i t ( R2r - r 3 / 3 ) / V = Nr(3R x - r 2)/4R 3 .

For a quadrant, Ve = ~ R 3 / 3 , and there are two planar faces, each of area 7rR2/2. The expected number of atoms within r of each of these planar faces is given by E[ Nt, l = N r ( 3 R 2 - r 2 ) / 8 R 3 .

gP ~ = 7rZ(2r3/3

=

V-----s

FI

If we assume that atoms are uniformly distributed within a cluster, their distances from the planar faces will also be uniformly distributed. Thus, if there are N e atoms within r of the planar face, men the expected distance from the planar face of the ith atom is di = i r / ( Nt, + l ) ,

3

/ 0

1r( R 2 - x 2 ) d x

An octant has volume zrR3/6 and three planar faces, each of area ~rR2/4. The expected number of atoms within r of each of these planar faces is given by E[ Np] = N r ( 3 R 2 - r 2 ) / 1 6 R 3.

Using these expressions for E[ N/, ] and assuming an equal probability of applying each of the three variants of the X2 crossover operator, Fig. 8 shows how E[ Arc] varies as a function of R (N = 40, K = 20). Clearly the number of close atoms generatexi is considerably less than that generated by the X! crossover operator (Fig. 7). 4.3. Analytical-experimental comparison

As a means of confirming the expressions obtained above for E[ Nr ] and E[ Np ], optimisation of a 40 atom cluster was performed and the tbllowing data recorded: • the actual number of close atoms after ,~'~ and X2 crossover operations;

146

W.J. Pulkm/ Computer Physics Communications 107 (1997) 137-148 ,

,

- -

,

_

T

,

11

16

111

0.3

A

I~ 11

14

g

N

II

~ 0.25

M

E ,,

Expecl(~.l

Average

&a

&

]

~0.15

i°

~

0.1

W

O.O5

o

IN

o

II

II

II

II

"5

5 Generation

(a) X1 - expected and average close atoms.

Ill

,'o

Generation

,'s

(b) X2 - expected a n d average close atoms.

Fig. 9. Expected and average number of close atoms for XI and 2'2 crossovers. The expected number is calculated by determining the average radius of the clusters in the pool and using the expressions derived above for El Nrl and El Npl. The observed average number of close atoms is determined directly from the child clusters.

® the value of R for each parent cluster. Averaging the close atom counts across child clusters for each generation produced the solid lines shown in Figs. 9a (Xt) and 9b (X2). Using the average values obtained in a single run for R, in conjunction with the expressions obtained above for E[ Nr ] and E[ Nc ], produced the data points flagged by stars in these plots. The sudden increase in the average and expected number of atoms was caused by a significant reduction in R for the population between generations six and seven. Given the stochastic nature of the process and the approximations used in the analysis, there is reasonable correlation between theory and actual results. Clearly the 2'2 geometric crossover operator is particularly effective at generating new clusters with a minimal number of close atoms. As shown by Fig. 6, this leads to a greater probability of a lower final energy in the subsequent BFGS optimisation. 4.4. X2 crossover enhancement

From the results presented, crossover operators which tend to produce lower numbers of close atoms in child clusters are, on average, more effective than those that produce higher numbers of close atoms. For the X2 crossover operator, if there are no close atoms in the parent clusters, then the only possibility

for close atoms in the child clusters is at the boundaries of the parent cluster volume elements used to create the child clusters. To eliminate this possibility, the X2 crossover operator can be modified (X~) in the following way to guarantee that it produces child clusters with no close atoms: For the hemisphere variant of X2, opposing hemispheres (above and below the xy-plane) from each parent cluster were initially translated directly away from the xy-plane and then moved together, in small incremental steps, with the cluster energy calculated at each step. Child clusters constructed from quadrant and octant volume elements were created using an analogous procedure. As a proportion of local minima occur at high energies [ 24], this physical translation of atoms reduces the complexity of the optimisation problem by effectively eliminating these local minima. In fact, Neumaier [24] estimates that removing these high energy local minima reduces the number of local minima on the PES from O(exp ( N 2 ) ) to O(exp (N) ). Other global optimisation methods for atomic clusters have eliminated these high energy local minima using either a physical separation of closc atoms [ 2 ] or a uniform expansion about the centre of mass of the complete cluster [ 25 ].

W.J. Pullan/Computer Physics Communications 107 (1997) 137-148

oI

,

,

,

,

,

•

,

,1"

,

',

147

,

,

,

~ o

-20

g

-2O

.m

../-

.E_ - 4 0

--4O o

r,f}

o

~-~

L-so

;r!

~-~

: . : / f •

.,~ - 1 0 0

8

-120

/ -1~14t

;

." ,.-..

-120

,/ i -120

I -100

- 8 0' - 6 0' - 20 Average patenl cluster energy

- 2'o

0

(a) X'~ - energy after crossover.

-140 -141

- 1 2'0

- 1'0 0

- ~0

Averageparent

'

-60 cluster energy

20

-

2'o

0

(b) X~ - energy after BFGS optimisation.

Fig. I 0. Average parent, crossover and BFGS optimi~d cluster energies for the X~ crossover operators during the optimisation of an N = 30 atomic cluster. All positive energy clusters have been omitted from this plot. Each point represents a single BFGS optimisation and shows the relationship between the average parent cluster energy and either the child energy after crossover or the final optimised energy of the child cluster. Points below the solid diagonal line show where the child cluster had a lower energy than the average parent cluster energy.

Fig. 10 shows the average parent, crossover and BFGS optimised cluster energies for the X~ crossover operators during the optimisation of an N --- 30 atomic cluster. For X~, 7% of child clusters had an energy lower than the average parent cluster energy immediately after crossover while 32% had an energy lower than the average parent cluster energy after a subsequent BFGS optimisation (as compared to 5% and 21% respectively for X2 and 0% and 13% for Xl ).

5.

Conclusion

The GA described in this paper was successful in solving atomic cluster problems ( N = 2 . . . . . 88) for the followir~g reasons: ® in common with all GAs, it explores the search domain in parallel using all atomic clusters in the current population as starting points; ® the use of the BFGS local optimisation to reduce the search space to local minima only; ® the physical separation of volume elements during the X~. crossover operation to ensure that there are no close atoms in the child clusters (which effectively eliminates all high-energy minima from the search domain);

• the GA operates at a number of levels when generating new atomic clusters. At the most coarse level is the hemisphere version of g~ which interchange large volumes between current clusters. Fine tuning of the cluster is provided by the quadrant and octant versions of g~ and the use of BFGS and APSE optimisation within the Mutate operator. The major constraint on the future development of this GA is the processor requirements of the BFGS local optimisations. Improvements in performance will come from a better balance of processor utilisation, implementation of more effective methods for early termination of the BFGS local optimisation and refinements to the crossover operators. Investigations are currently being performed so that these improvements may be implemented.

Acknowledgements

The author would like to thank Graham Wood for his many critical and constructive comments in preparing this paper. For a substantial portion of this study the author was a visitor at the Computer Science Department, University of Colorado at Boulder, and would particularly like to thank Robert Schnabel, Richard Byrd and Elizabeth Eskow for the opportunity to work

148

W.£ PullanlComputer Physics Communications 107 (1997) 137-148

with their g r o u p and for sharing their k n o w l e d g e o f global optimisation o f molecular structures.

References l l l D.M. Deaven and K.M. Ho, Phys. Rev. Lett. 75 (1995) 2. [21 D.M. Deaven, N. Tit, J.R. Morris and K.M. Ho, Phys. Chem. Lett. 256 (1996) 195. [3] W.J. Pullan, J. Comput. Chem. 18 (1997) 1096. [4] WJ. Puilan, Ph.D. Thesis, Department of Mathematics and Computing, Central Queensland University (1996). J5] M.R. Hoare, Adv. Chem. Phys. 40 (1979) 49-135. [6] M.R. Hoare mad J.A. Mclnnes, Fara. Disc. Chem. Soc. 61 (1976) 12-24. [7] M.R. Hoare and P. Pal, Adv. Chem. Phys. 20 ( 1971 ) 161196. 181 M.R. Hoare, J.A. Mclnnes, Adv. Phys. 32 (1983) 791-821. 191 J.A. Northby, J. Chem. Phys 87 (1987) 6166-6177. !1Ol G. Xue, J. Global Opt. 4 (1994) 2425-2440. [ I l l L.T. Wille, Chem. Phys. Lett. 133 (1987) 5. 1121 G. Xue, J. Global Opt. 4 (1994) 2187-2208. ll31 C.D. Maranas and C.A. Floudas, J. Chem. Phys. 97 (1992) 7667.

ll4l J. Kostrowicki, L. Piela, B.J. Cherayii and H.A. Scheraga, J. Phys. Chem. 95 ( 1991 ) 4113-4119.

1151 R.J. Wawak, M.M Wimmer and H.A. Scheraga, J. Phys~ Chem. 96 (1992) 5138-5!45. ll61 B. Hartke, J. Phys. Chem. 97 (1993) 9973-9976. 1171 J. Mestres and G.E. Scuseria, J. Comput. Chem. 16 (1995) 6. 1181 S.K. Gregurick, M.H. Alexander and B. Hard~e, J. Chem. Phys. 104 (!996) 2684-2691. I191 D. Goldberg, Genetic Algorithms in Search, Optimization and Machine Learning (Addison-Wesley, Reading, MA, 1989). 1201 J. Antonisse, Proc. Third Int. Conf. Genetic Algorithms, George Mason University (Morgan Kaufmann, Los Altos, CA, 1989). 121l C.Z. Janikow and Z. Michalewicz, Proc. Fourth Conf. Genetic Algorithms (Morgan Kaufman, San Mateo, CA, 1991). 1221 K.P. Wong and Y.W. Wong, Proc. ANZIIS-93 (Perth, Western Australia, 1993). 1231 W.J. Puilan, Aust. Comp. J. 29 (1996) 4. 1241 A. Neumaier, SIAM Rev., to appear. 125l R.H. Byrd, T. Derby, E. Eskow, K.P.B Oldenkamp and R.B. Sehnabel, CU-CS-652-93 (Dept. of Computer Science, Univ. of Colorado, Boulder, i 993).