Unconstrained optimisation in surface crystallography by leed: Preliminary results of its application to CdTe(110)

Surface Science 187 (1987) 175-193 North-Holland, Amsterdam

175

U N C O N S T R A I N E D O P T I M I S A T I O N IN S U R F A C E C R Y S T A L L O G R A P H Y BY LEED: P R E L I M I N A R Y R E S U L T S OF ITS APPLICATION TO CdTe(ll0) P.G. C O W E L L Department of Physics, University of York, Heslington, York Y01 5DD, UK

and V.E. D E C A R V A L H O Departamento de Fisica - ICEx, Universidade Federal de Minas Gerais, Belo Horizonte, Brazil

Received 21 November 1986; accepted for publication 15 April 1987

This paper describes the application to LEED of a form of the well-known method of unconstrained optimisation, originallydeveloped by Hooke and Jeeves. It has been found to be of great use in the iterative searching that must be performed as a part of the process of discovering the atomic geometry of a surface using LEED. By using such an optimisation scheme, the speed with which a crystal surface structure may be found is greatly increased. The codification of the search strategy used in LEED in this way should also result in greater confidence being placed in LEED as a surface structure analysis technique. The CdTe(ll0) surface is used as the first application of this new scheme using newly collected, high quality experimental data. Nine structural and non-structural parameters were investigated to give a preliminary result for this surface.

1. Introduction The purpose of this paper is to propose a new procedure for use in surface structure studies by LEED, which may possibly result in an e n h a n c e m e n t in the reliabifity of structures determined by this method, b u t certainly improves the speed and efficiency with which a result may be obtained. This new method involves a n application of what is, by present day standards, a rather unsophisticated form of u n c o n s t r a i n e d optimisation to the modelling stage of a L E E D investigation. The advantages in terms of efficiency a n d time saving will be made evident in general terms. The possible advantages to be gained from the systematisation of search strategies by making the exploration of a parameter space more rigorous will be discussed. 0 0 3 9 - 6 0 2 8 / 8 7 / $ 0 3 . 5 0 9 Elsevier Science Publishers B.V. ( N o r t h - H o l l a n d Physics Publishing Division)

176

P.G. Cowell, V.E. de Carvalho / Unconstrained optimisation in crystallography

This paper is concerned with the theoretical or computer modelling part of LEED. Recent improvements made in LEED diffractometers have significantly decreased the time taken by the experimental phase of LEED. Consequently, efforts are now being made to speed up LEED calculations mostly by what are often called "approximate" methods. These include techniques such as quasi-dynamical LEED [2,3], tensor LEED [4], beam set neglect and others [5,6]. This paper, however, describes an increase in efficiency which allows greater speed in the theoretical part of LEED without needing to perform the individual calculations more quickly. The theoretical part of LEED consists both of the calculation of LEED data, (I(V) curves) and the comparison between theory and experiment, usually using a dissimilarity criterion (an R-factor) and the construction of the theoretical model to be used. By using an optimisation technique, it should be possible to reduce the number of calculations needed to solve a structure by intelligent use of the results of the previous comparisons. The Hooke and Jeeves method of unconstrained optimisation [1] has the properties required for use in LEED. More sophisticated optimisation schemes [7] which need prior knowledge of the shape of the space being explored are inappropriate in LEED where such knowledge is unobtainable. A recent LEED controversy arising from the way in which a solution is searched for has involved the GaAs(ll0) surface. The structure was solved using LEED in 1976 [8] but an alternative structure proposed by isochromat spectroscopy [9,10] prompted a re-evaluation of the earlier LEED result [11]. A possible new solution was found as demonstrated by a previously undiscovered minimum in the R-factor. Evidence has been produced to show that this is not the best solution and, in fact, the original surface structure gives a minimum in the R-factor deeper than the newer proposed structure. This has generated some discussion in the form of papers criticising the way in which the surface structure should be derived from experimental data [12-14]. The work done on GaAs(ll0) has been generally of high quality in terms of the experimental data, the theoretical calculations and the amount of data used in the comparisons between theory and experiment. The controversy rests almost entirely in the way in which the complex parameter space was searched, i.e. the range of trial theoretical models used. Similar results and potentially similar controversies are provided by nearly all of the (110) surfaces of the zinc-blende compound semiconductors which have been found to have very similar atomic geometries [15]. If a systematic search of a wide part of the parameter space had been made in the early investigation of GaAs(ll0), the recent disagreements over this surface might have been avoided. Although this may be said with hindsight, it also should be pointed out that it is only recently that it has been thought reasonable to invest computer time in exploring such a complicated surface in great detail [14]. This paper will now describe a 'search strategy for use in

P.G. Cowell, V.E. de Carvalho / Unconstrained optimisation in crystallography

177

LEED surface modelling based on an application of the Hooke and Jeeves unconstrained optimisation scheme.

2. Description of the search strategy The way in which a LEED investigation- of a surface is performed may be described as a search to minimise the differences between theoretical and experimental data, as measured by an R-factor, through varying the model used to calculate the theoretical data. The search is carried out in a multidimensional parameter space consisting of a mapping of an R-factor on the values of the structural and non-structural parameters which define a crystal surface. The aim of the search is to find a minimum efficiently, i.e. using the minimum number of surface model calculations. The Hooke and Jeeves type of search method [1] uses information about the shape of the R-factor function in the immediate vicinity of a point in the parameter space (called a base point) to determine the best direction in which to move in order to reduce the value of the R-factor function. The search is then moved in this direction until the value of the function ceases to diminish. When this point is reached, a new base point is set up and the shape of the function is assessed again to find the new direction in which the function decreases most rapidly. New base points are repeatedly set up in tiffs way moving closer and closer to the position of the minimum. This method of searching is more efficient than the simple search method of minimising each parameter value individually and successively. By changing several parameter values simultaneously, the R-factor function may be made to reduce more quickly than by varying any parameter in isolation. Also it is well suited to the way in which LEED structures are investigated. The Hooke and Jeeves search technique consists of two distinct parts: the investigation of the area around each base point and the moves made from the base point along the best direction for reducing the R-factor values. The way in which these two operations have been implemented for use in L E E D will now be described.

2.1. Exploration For each of the parameters which are to be varied in the search there is an associated step size. In the present work, these step sizes were chosen by physical intuition and the experience of investigations of similar surfaces. Two values for each parameter are considered in the exploration, one step either side of the base point parameter value. The R-factor for each parameter is then minimised and all the new parameter values combined into a new model,

178

P.G. Cowell, V.E. de Carvalho / Unconstrained optimisation in crystallography

the R-factor for which is expected to be lower than that at the base point or any of the exploratory points. 2.2. Move

The search keeps moving in the direction deduced from the exploration in equal steps until the R-factor fails to decrease. At this point a new base point is established and a new exploration initiated. Eventually it should be found that no move is possible. The exploration indicates a minimum R-factor at the base point. A decision is then needed: the step sizes for the parameters may be reduced or, if it is considered that they are already sufficiently small, the search may be considered complete. A certain amount of care is necessary when choosing the initial base point. It is clear from the way in which the search always proceeds downhill from the base point that, in a space with several minima, the minimum chosen as the final goal of the searching technique will be the one nearest to the initial base point. It is essential when using the method as a tool in a LEED investigation that some information on the shape of the R-factor function should be obtained before choosing the starting point for the search. This information may come from earlier LEED investigations which it is hoped may be refined using better experimental data or computational techniques. Or models proposed by surface structure techniques other than LEED may be used. The best source of information is provided by a rough survey of the R-factor function achieved by varying a limited number of parameters (those expected to be the most significant) in large steps spanning the whole range of parameter space where the solution could reasonably be expected to be. This expects judgements to be made about which are the most significant parameters and what limits may be put on their possible values. This may not be possible in all cases and these a-priori decisions should always be confirmed by observing the way in which a search proceeds as they may effect the validity of the final result. No serious attempt to address the considerable problems associated with the bounding of the parameter space in a LEED search will be made in this paper. It is a problem that is very difficult to generalise. General physical principles may be used as guides to the limiting values of parameters; knowledge about the bonding between atoms can give an idea of the maximum and minimum distances allowable in a crystal; theoretical solid state and atomic physics can provide information about the non-structural parameters in a LEED calculation. Probably the most useful and potentially the most misleading source of information is previous experience of LEED investigations, particularly of similar types of surfaces. This is useful because it may be used to speed up an investigation considerably, but'it may prove misleading if,

P.G. Cowell, V.E. de Carvalho / Unconstrained optimisation in crystallography

179

by taking short-cuts, the investigation is prejudiced in favour of a structure that fits in with previous results. These problems of bounding the parameters in the search and of choosing the significant parameters apply only the initial survey. The parameters in the Hooke and Jeeves type of search are unconstrained. The minimum immediately downhill from the initial base point will be found if one exists. It is mostly in the interpretation of the result of the search that the information

I n i t i a l s u r v e y to f i n d suitable starting

points

for

the

the

best

search.

Explore

to

find

Reduce the step s i z e s used in the e x p l o r a t i o n .

d i r e c t i o n in which to move.

, Is any

move

possible?

I Move in until

the

the

'best

R-factor

I

~_~

I

Is this solution?

an

acceptable

direction' ceases

to

reduce.

Fig. 1. Simple flow d i a g r a m o f the H o o k e and Jeeves t y p e o f optimisation scheme described in this paper.

180

P.G. Cowell, K E. de Carvalho / Unconstrained optimisation in crystallography

gleaned from the initial survey is important. Unless it can be shown that the minimum or minima found be a search are the only important ones in the investigation, there will remain a lingering doubt over the result. Fig. 1 shows a summary flow diagram of the proposed search method.

3. Application to CdTe(ll0) The CdTe(110) surface was chosen as a test case for the development of the Hooke and Jeeves method as a LEED searching technique. It is a complex surface in which the top two layers are expected to reconstruct away from the bulk structure with different relaxations for the two atom types. The symmetry of the LEED pattern indicates that the (110) surface retains the one mirror plane of the bulk (110) planes. An extensive investigation of the zinc-blende compound semiconductors including CdTe has been made in the past [15,16] that provides information valuable when evaluating this new technique. This area of study, the zinc-blende semiconductors, has not been without controversy [12,13] as indicated in the introduction to this work. LEED alone has not been particularly successful in distinguishing between the two minima in R-factor found in the parameter space explored for GaAs(ll0) [17]. The less-favoured minimum is reported not to be as deep as the other but this subsidiary minimum has rarely been investigated in detail. Thus the CdTe(ll0) surface offers the advantages of a moderately complex structure which has great similarities to structures already explored. Something is already known about the shape of the R-factor as a function of the parameters from the earlier investigations of CdTe(110) and similar surfaces. In particular two minima are expected [14]. This provides an interesting test for the search technique. Also a more detailed exploration of the subsidiary minimum than has been attempted previously is now possible without inordinate amounts of time and effort.

3.1. Experimental The experimental data was collected using a computer controlled electron diffractometer [18]. The details of the collection of the data very similar to those for work already published on InSb(110) [19]. A cubic single crystal of CdTe (about 1 • 1 x 1 cm 3) was cleaved in vacuo (5 • 10 -11 mbar) at room temperature. A flat area about 0.5 cm 2 of (110) surface was produced. This area was examined using an in-situ UHV scanning electron microscope to probe the surface topology and Auger analysis was used to examine the surface chemically. It was found to be flat and uncontaminated. A clear diffraction pattern with the expected symmetry was obtained. A great amount of effort was put into determining experimentally the angle

P.G. Cowell, V.E. de Carvalho / Unconstrained optimisation in crystallography

181

of incidence of the electrons on the surface. Because the diffraction pattern displays symmetry with only one mirror plane, the standard technique of adjusting the angle of incidence until the symmetrically equivalent beams are degenerate is not sufficient in this case completely to define normal incidence, except in a single plane perpendicular to the surface. The method used to set up the best approximation to normal incidence has been described elsewhere [19] as it was applied to InSb(110). Twenty beams were collected nearest to normal incidence. They were all corrected for variations in the current of incident electrons. Degenerate beams were averaged together and twelve of the resulting beams were used for comparison with the theory. The experimental I ( V ) curves were smoothed once using a nine-point averaging algorithm. Experience performing the experiment led to doubt being cast on the validity of data collected below 50 eV. No use was made of this experimental data when calculating R-factors. Only one angle of incidence, 0, was varied when searching for a solution. It is a reasonable assumption that the angle of incidence will be close to zero in the mirror plane of the sample (0) and almost exactly zero out of the mirror plane. 3.2. Calculational details

All the calculations for CdTe(ll0) described in the following sections were performed using a symmetrised version of the CAVLEED program [20] capable of multiple scattering calculations involving two atoms per unit cell. The crystal potential was a muffin-tin potential calculated from the superposition of spherically averaged wavefunctions taken from the calculations of Clementi and Roetti [21]. Phase shifts were calculated for nine angular momentum values at fifteen energies up to 10 hartrees. The exchange parameter a = 0.6667. 59 symmetrically inequivalent beams were used in the intralayer calculation of transmission and reflection matrices. Renormalised Forward Scattering [22] was used for the interlayer scattering calculation. When this failed to converge, as occasionally happened, the program switched to the layer doubling method. The beam intensities were computed for incident energies from 50 to 120 eV in steps of 5 eV. The calculations were performed on the CRAY-1S at the University of London Computing Centre. 3.3. CdTe(l lO) - initial survey

The importance of some knowledge of the way in which the R-factor varies with parameters has already been pointed out. The search technique may almost be regarded as a technique for refining the information from the initial survey. The result of the search is to find the minimum nearest to the starting point so some idea of the probable result, however crude, most be obtained if

182

P.G. Cowell, V.E. de Carvalho / Unconstrained optimisation in crystallography To

Cd

Fig. 2. View of the zinc-blende (110) surface from the (~10) direction. On the left of the picture the bulk structure is shown and on the right a reconstructed surface. The parameters varied in the searches are shown. They are measured with respect to the bulk structure in the directions indicated by the arrows. A negative value for a rumple (either 01 or 02) implies greater movement of the tellurium atoms towards the vacuum than the cadmium atoms.

a reliable result is to be obtained. Experience with other zinc-blende (110) surfaces indicates that ~0] and q (defined in fig. 2) will be most significant in determining the structure. A contour plot of the Pendry R-factor as a function of these parameters is shown in fig. 3. There are clearly two minima. One is at a rather large bond rotation angle, ~01, and this corresponds closely to the structure published by Duke et al. [23]. The other minimum is shallower and corresponds to a bulk termination, ~1 = 0 and e a = 0. Previous work and the initial rough survey establish that there are two minima to be considered indicating the need for two searches. They also give two reasonable starting points for the searches. One at the model published by Duke et al. [23] which is expected to prove to be the correct solution and the other at a bulk termination. 3.4, Details of the search Eight parameters were varied in each of the searches: top layer rumple p], top layer relaxation q , top layer anion shift Yal' top layer cation shift Yc,, second layer rumple P2, second layer relaxation ~2, angle of incidence O, and Debye temperature T. The meanings of these parameters are explained in fig. 2. In this paper all the relaxations, rumples and the lateral shifts of the cations and anions will be given as percentages of the bulk lattice constant, 6.482 A. Angles will be in degrees and temperatures in kelvins unless stated otherwise. This set of parameters was chosen principally from the experience of previous work [14,19]. The Debye temperature was the same for all the atoms in the calculations. There is some evidence [24] that, in the surface of zi'nc-blende semiconductors,

P.G. Cowell, V.E. de Carvalho / Unconstrained optimisation in crystallography

l

183

~ . - l - -

-5.0

.~...:U-__1o ,~,/o

~

~,/~>

5.o

-20

-,o

b

2_S

~

_I2

~'~

(/~

k

'

;

"i

',/%

0.8~'~-20

:,/j '

/

1o\ /

~)

98

" -4o"

0)lJ/o

Fig. 3. (a) Quasi-three-dimensional plot of a function of R-factor against the top layer relaxation, cl, and the top layer bond rotation angle, ~1. The meanings of these symbols are given in fig. 2. The function, F, of the R-factor, R, is: F(R) = (R - R min + 0.1)-2, where R mi, is the minimum value of R in the range of values plotted. This function of R serves to emphasise any minima in the R-factor surface by transforming them into peaks that are much sharper than they would otherwise be. In this figure two minima are indicated, the deeper of them being at the edge of the plot ( q = - 5.0%). (b) Contour plot of the R-factor data presented in (a) indicating the actual values of the R-factors for the two parameters varied (( t and Wl). The two minima are at c 1 = - 5.0%, wt = - 3 0 ~ (the deeper minimum) and c1 = 0.0, ~1 = 0.0.

the atoms have bulk-like vibrational amplitudes. The Debye temperature was i n t r o d u c e d i n t o t h e c a l c u l a t i o n s t h r o u g h c o m p l e x p h a s e shifts. U p to n i n e phase shifts were used for each energy point. T h e R - f a c t o r c h o s e n f o r u s e i n t h e s e a r c h w a s t h e P e n d r y R - f a c t o r [25]. T h i s R - f a c t o r h e a v i l y w e i g h t s t h e p o s i t i o n s of t h e p e a k s a n d t r o u g h s i n t h e I ( V ) curves when assessing their similarity. The relative peak heights are not

184

P.G. Cowell, V.E. de Carvalho / Unconstrained optimisation in crystallography

considered of major importance in this measurement. Although the Pendry R-factor was used as the similarity criterion in the search, two other R-factors were calculated to provide extra information on the way in which relative peaks heights changed. They were the Euclidean N o r m and the Weak Integrated Distance [26] both of which are sensitive to the relative peak heights or envelope of the curves being compared. Both are metric distances (see ref. [27]). The Euclidean N o r m is the integral of the root mean square difference between two curves both normalised to unit integral over the range of comparison. A function is a norm if, in addition to being a metric, it is linear. The Weak Integrated Distance is defined as follows:

R(f, g)= f o l l F ( x ) - G ( x ) l dx, where F and G are nondecreasing functions of the from zero to one:

F(x)=

So;

( t ) dt

and

G(x)=

I(V)

curves f and g rising

/0xg(t) dt.

Both f(t) and g(t), the experimental a n d / o r theoretical normalised to unit integral.

I(V)

curves, are

3.5. CdTe(110) - the search The first search made was to find the position of the main minimum (the one corresponding to a surface bond rotation near 30 o ). The initial point for the search was the structure published by Duke et al. [23]. The structure is defined in terms of the parameters used in this work in column 1 of table 1. The errors in the table are not those published by Duke et al. [23] but are Table 1 Structures from first search

Pi (%) ~1 (%) Ya~ (%) Yc, (%) 02 (%) c 2 (%) 0 (deg) T(K) R-factors (shift (eV)) Pendry Euclidean Weak distance

Initial model (column 1)

Final model (column 2)

-12.7• -3.6• 6.0• 9.4 2.8 0.0_+1.7 0.0• oo

-12.4• 1.3 - 3 . 6 + 1.4 5.2+ 1.9 8.6• 3.8 2.5_+ 3.6 0.5• 1.1 -0.5• 1.7 350 +300

0.5640 (1) 0.1780 (3) 0.0937 (3)

0.5189 (0) 0.1741 (1) 0.0853 (2)

P.G. Cowell, V.E. de Carvalho / Unconstrained optimisation in crystallography

185

Table 2 Structures from second search Initial model (column 1)

Final model (column 2)

~1 (%) Ya, (%) Ycl (%) P2 (%) c 2 (%) 0 (deg) T (K)

0.0-t-2.1 0.0+2.0 0.0 0.0 + 4.9 0.0+3.1 0.0 + 1.9 0.0 + 3.1 oo

0.0+1.9 0.0+1.8 0.0 0.0 + 4.4 -2.0+2.0 0.0 + 1.7 2.0 + 2.3 oo

R-factors (shift (eV)) Pendry Euclidean Weak distance

0.7481 ( - 4) 0.2302 ( - 3) 0.1062 (1)

Pl (%)

0.6344 ( - 2) 0.2206 ( - 2) 0.1189 (1)

d e r i v e d f r o m the p r e s e n t w o r k as d e s c r i b e d b e l o w . I n all the c a l c u l a t i o n s t h e r e was o n e e x t r a p a r a m e t e r , the real p a r t o f the i n n e r p o t e n t i a l (the i m a g i n a r y p a r t was fixed at - 4 eV). T h i s p a r a m e t e r was n o t u s e d in t h e s e a r c h as t h e r e is n o c o n s t r a i n t o n the cost of its calculation. It is well a p p r o x i m a t e d b y a shift in the t h e o r e t i c a l curves. T h i s shift was e v a l u a t e d for e a c h t h e o r e t i c a l m o d e l i n d i v i d u a l l y as t h e R - f a c t o r was m i n i m i s e d . T h e R - f a c t o r s for the initial m o d e l are also s h o w n in c o l u m n 1 of t a b l e 1 t o g e t h e r w i t h the shifts w h i c h g a v e the m i n i m u m R - f a c t o r values. A f t e r f o u r i t e r a t i o n s o f the s e a r c h l o o p t h e r e c e a s e d to b e a n y s i g n i f i c a n t i m p r o v e m e n t in the P e n d r y R - f a c t o r e v e n t h o u g h t h e sizes of the steps w e r e

2

-

0.80

1

o

..__....~-'~

'

_t-~ . . . . . . .

-1

-2

-3

~

,

.-

-r

. . . . . . . . . .

2.0

3.0 O/~

-

_

; ...... ~176

~ "~

-4

9..

9.

,., ,., ~n.~=

,...~

.~

,

. 9

~176176

Fig. 4. Contour plot of Pendry R-factor against P2 (see fig.2) and 0. The minimum is at P2

=

- -

2.0%, 0 = 2.0 o.

186

P.G. Cowell, V.E. de Carvalho / Unconstrained optimisation in crystallography

reduced to quite small values, smaller in fact than the error in the parameters that could be derived from the R-factor. The final result of the search, the best model, is summarised in column 2 of table 1 together with the corresponding R-factor values. The subsidiary minimum near a bulk termination was then investigated starting from the bulk values. All the parameter values were zero except the Debye temperature which was infinity, indicating a rigid lattice. (See table 2, column 1.) I(V)

arbitrary

units

Experiment

v 0

!

I

100

200

eV

Fig. 5. Experimental and three calculated curves for the (11) beam. This has a low individual beam Pendry R-factor of 0.24 (curve 2). Theory curve 1 is calculated for the model of Duke et al. [23] (the initial model for the first search). Curve 2 is the final model for the first search, the best result (see table 1). Curve 3 is the final model for the second search (table 2).

P.G. Cowell, V.E. de Carvalho / Unconstrained optimisation in crystallography

187

The search again converged quickly on a minimum: in four iterations of the technique. The result of this second search is shown in column 2 of table 2. Only two of the parameters have changed from their bulk values, P2 and 0. A contour map of R-factor plotted against P2 and 0 in fig. 4 shows the minimum. An example showing the experimental and various calculated I ( V ) curves may be seen in fig. 5.

4. Errors

The search technique described here is well suited to the calculation of errors in the parameters from the Pendry R-factor values. The evaluation of the error associated with a parameter depends on a knowledge of how the R-factor varies as the parameter is changed about the minimum value. The R-factor calculation can provide something analogous to an error in R-factor via the calculation of RR-factors (see ref. [25]). Only parameter values for which the difference between its associated R-factor and the minimum R-factor is less than the R-factor "error" are considered to be significant. At the end of an exploratory stage of the search, the R-factor data is ideal for use in calculating errors. A parabola may be fitted through the three R-factor values for each independently varied parameter and used to give the range of error for each parameter. This type of calculation is described pictorially in fig. 6.

R-factor

t'-'error"in R-factor

I

1-,--

I

"-~[

,

Parameter

r a n g e of e r r o r in p a r a m e t e r

Fig. 6. Diagram showing schematically the way in which the error in a parameter is evaluated from a Pendry R-factor calculation.

188

P.G. Cowell, V.E. de Carvalho / Unconstrained optimisation in crystallography

All the errors quoted above were calculated in this way from the Pendry R-factor. Errors may only be calculated if a parabola with a minimum may be fitted through the three R-factor values. When far from the optimum solution, errors could not always be calculated. It was found that the effect of reducing the step size even to values substantially lower than the error found for the parameters had little effect on the resultant errors. The correct parabola could be fitted through the three R-factors even when the corresponding parameter values were very close together. This indicates that the sensitivity of an R-factor to changes in a parameter is very much greater than the errors quoted on that parameter may lead one to believe [28].

5. Discussion of CdTe(ll0) structure This paper has just presented a new provisional structural determination of the CdTe(ll0) surface. However, the main purpose of this paper is rather to describe the technique used to find this solution. In this section the significance of the result as a surface structure will be addressed in brief. The result presented can confidently be taken as a reasonable representation of the structure of CdTe(ll0). It confirms the result of Duke et al. [23] using independently collected experimental data and different multiple scattering calculations. The slight differences between the present model and that of Duke are consistent with the analogous results found for InSb(ll0) [19] in recent work. InSb and CdTe are very similar (isoelectronic) substances and, despite differences in the ionicity and therefore the nature of the chemical bonding, their (110) surface structures have been found to be very similar [23] and to make sense in terms of bond rehybridisation and energy minimisation

[24]. The presence of a subsidiary minimum in searches of zinc-blende (110) surfaces has been regularly noted [11,14,17,29]. In this work, for the first time for CdTe(ll0), the immediate area in the parameter space around this subsidiary minimum has been carefully explored. It is possible to say with confidence that, in the space defined by the parameters used in the search, the subsidiary minimum is shallower than the main minimum found for higher surface bond rotations. The difficulty of assessing the relative depths of minima by varying only one parameter (i.e. by taking cuts through the parameter space) has led to misleading results in the past [14]. Even though the structure as described here is a reasonable structural result, it is not expected ultimately to be the best obtainable. The work is not yet complete. There is more experimental data collected away from normal incidence from which it is yet hoped to elicit greater precision in a structural result. The Debye temperature chosen by the Pendry R-factor was signifi-

P.G. Cowell, V.E. de Carvalho / Unconstrainedoptimisation in crystallography

189

cantly and systematically different from that which would have been chosen by either of the other two R-factors. The Pendry R-factor consistently chose very high Debye temperatures whereas the others chose much lower values nearer the bulk Debye temperature. The variation of the R-factors over a small range of Debye temperatures may be seen in fig. 7. The reason for this discrepancy is that while high Debye temperatures result in I ( V ) curves where the positions of the peaks are accurately reproduced the relative peak heights are quite wrong. This results in a low Pendry R-factor but high Euclidean N o r m and Weak Integrated Distance values. The relative peak heights are well reproduced if a low Debye temperature is used. The Euclidean N o r m and Weak Distance values are therefore much lower but the peak positions in the I ( V ) curves have moved relative to the experiment and the Pendry R-factor is consequently high. It has been suggested that the surface may have a different and higher Debye temperature than the bulk. This possibility will be investigated.

R-factor 0.534

0.530

0,526

0.522

0.08,!

0.083

0.081 I

I

I

250

300

350

9

Debye temperature/K

Fig. 7. R-factor values plotted as a function of Debye temperature. The top curve, (a), shows the Pendry R-factor decreas!ngwith increasing Debye temperature. This is in conflict with curve (b) which shows that Weak Distance R-factors have the opposite trend. In both cases three values are plotted and the only parameter to change was the Debye temperature.

190

P.G. Cowell, V.E. de Carvalho / Unconstrained optimisation in crystallography

6. Discussion of search strategy The search was able to find two minima and to distinguish the better of them. In the case of CdTe(ll0) it is very unlikely that, if there are any other minima they will be anything but very shallow local minima. Pl and cl were found to be the most significant parameters in determining the surface structu(e, vindicating their choice for use in the initial survey. There are only two R-factor minima in the parameter space of these two parameters and only these are expected to be significant, i.e. possible solutions. In general, information about the relative sensitivity of the parameters will not be available a-priori. It must be evaluated during a search through parameter space as each of the parameters is varied in turn. The Hooke and Jeeves search is a method of unconstrained optimisation so (during the search) there is no constraint placed on the values that the parameters may take. This is not true of the initial survey stage. It must be performed over a finite space but the problem of deciding the range of that space cannot in general be solved easily. In the present study of CdTe(ll0) all the experience of the past has been used to determine the limits of the initial survey and although there are small signs of more minima to be found beyond the edges of the survey, they have not been considered important as the models to which they correspond are physically quite unreasonable.

6.1. Timings One of the main reason for developing a Hooke and Jeeves type search for use in LEED is to increase the speed with which a surface structure determination may be achieved. For the particular case of CdTe(ll0) presented here, there is the opportunity of comparing the time taken to solve two almost equivalent surface structures, one using the new technique and the other using a more intuitive search strategy. The other surface is InSb(ll0) [19]. In both cases the calculations were carried out using the same computing resources. All the LEED calculations were done on the CRAY-1S computer at the University of London Computing Centre and the calculation of R-factors, plotting of results, etc. was carried out using the DEC-10 computer at the University of York. The time taken to solve a surface structure by LEED may be measured in two different ways: either the total time for a conclusion to be reached or the computer time used in the process. The use of the new method of searching has shown improvements in the speed with which a structure determination is reached measured by either criterion. The comparison of total times is not easy to quantify. However, the new method, providing a logical scheme for the parameters to be varied systematically, produced a time saving of about a factor of five. Another measure of the time taken tb solve a structure is the

P.G. Cowell, V.E. de Carvalho / Unconstrained optimisation in crystallography

191

amount of CRAY computer resources used, as expressed in CPU executing time. For CdTe a total of 260 min were used, whereas InSb needed approximately 600 min. This includes the investigation of both potential structures for the case of CdTe whereas very little work was done for InSb looking at the subsidiary minimum in R-factor.

6.2. Possible improvements There is, therefore, a demonstrable increase in efficiency of a factor of at least three and possibly much more when comparing the new search strategy with the older intuitive way of searching. But several ways in which the search could be modified to make it more streamlined have suggested themselves which, together with speedier multiple scattering calculations, could enable far more complex structures to be investigated by LEED than could presently be contemplated. One way in which the search technique might be improved and made faster would be to interpolate between the three R-factors calculated for each parameter after the exploratory phase of the search. This has been attempted a-posteriori by fitting a quadratic function to R-factor values obtained early in the search. In many cases the minimum found in the quadratic was very close to the final parameter value. However, if the R-factor is insensitive to the variations in a particular parameter or the parameter values are far from the minimum R-factor the interpolation or extrapolation may be extremely inaccurate and misleading. As the search strategy is used now (as described in this paper), no use is made of the errors in the parameters calculated from the R-factor e~cept at the end of the search to indicate the accuracy with which the various parameters have been determined. It may, however, be possible to use either the errors calculated for each parameter or the errors in the R-factor calculated after each exploration in some way to determine future step sizes or to give a minimum significant change in R-factor.

7. Summary The search scheme described in this paper is an application of the well-tried Hooke and Jeeves method of unconstrained optimisation [1]. It may be summarised as consisting of three essentially separate processes. Firstly an initial survey is made of the shape of the R-factor over the parameter space to be used in the investigation. This gives the positions of the best starting positions for the search. The second and third processes in this scheme are an exploratory phase lOcated around a base point in the parameter space and a motion of the base point. The exploration is made in order to find the best

192

P.G. Cowell, V.E. de Carvalho / Unconstrained optimisation in crystallography

direction in which to move the base p o i n t so as to find the m i n i m u m m o s t efficiently. The searching m e t h o d is o u t l i n e d in the flow d i a g r a m shown in fig. 1. This scheme is suitable for d e v e l o p m e n t into an a u t o m a t i c searching technique because the j u d g e m e n t s a b o u t the way in which the theoretical model should be varied have been r e d u c e d to an algorithm. This will be of great value when used in c o n j u n c t i o n with the faster schemes being d e v e l o p e d for multiple scattering calculations [2-6]. Large areas of p a r a m e t e r spaces m a y then be explored very quickly. This systematisation of L E E D searches has the c a p a c i t y greatly to reduce the time taken to solve a surface structure, allowing for m o r e t h o r o u g h exploration of m o r e c o m p l i c a t e d p a r a m e t e r spaces to be c o n t e m p l a t e d . This can only serve to enhance the r e p u t a t i o n of L E E D as a technique for surface structure determination.

Acknowledgements The authors would like to t h a n k Professor M. P r u t t o n a n d Dr. S.P. T e a r for their helpful suggestions a n d discussions which have c o n t r i b u t e d to the development of the ideas p r e s e n t e d here. T h e w o r k d e s c r i b e d has been supported b y the Science and Engineering Research C o u n c i l ( U K ) . O n e of us (V.E. de C.) would like to t h a n k the C A P E S and U n i v e r s i d a d e F e d e r a l de Minas Gerais, Brazil for their financial assistance.

References [1] [2] [3] [4]

R. Hook and T.A. Jeeves, J. Assoc. Computing Machinery 8 (1961) 212. K. Heinz and G. Besold, Surface Sci. 125 (1983) 515. S.P. Tear, S.Y. Yousif and J.A.D. Matthew, Surface Sci. 175 (1986) L715. P.J. Rous, J.B. Pendry, D.K. Saldin, K. Heinz, K. Miiller and N. Bickel, Phys. Rev. Letters, submitted. [5] M.A. Van Hove and G.A. Somorjai, Surface Sci. 114 (1982) 171. [6] M.A. Van Hove, Approximations in LEED Intensity Calculations, Intern. Seminar on Surface Structure Determination by LEED and Other Methods, Erlangen, 1985. [7] For example, G.R. Walsh, Methods of Optimisation (Wiley, Chichester, 1985). [8] A.R. Lubinsky, C.B. Duke, B.W. Lee and P. Mark, Phys. Rev. Letters 36 (1976) 1058. [9] V. Dose, H.-J. Grossman and D. Straub, Phys. Rev. Letters 47 (1981) 608. [10] V. Dose, H.-J. Grossman and D. Straub, Surface Sci. 117 (1982) 387. [11] C.B. Duke, S.L. Richardson, A. Paton and A. Kahn, Surface Sci. 127 (1983) L135. [12] M.W. Puga, G. Xu and S.Y. Tong, Surface Sci. 164 (1985) L789. [13] C.B. Duke and A. Paton, Surface Sci. 164 (1985) L797. [14] P.G. Cowell, M. Prutton and S.P. Tear, Surface Sci. 177 (19,86) L915. [15] C.B. Duke, A. Paton and A. Kahn, J. Vacuum Sci. Technol. A1 (1983) 672.

P.G. Cowell, V.E. de Carvalho / Unconstrained optimisation in crystallography

193

[16] C.B. Duke, Atomic Geometry and Electronic Structure of Tetrahedrally-Coordinated Compound Semiconductor Interfaces, in: Surface Properties of Electronic Materials, Eds. D.A. King and D.P. Woodruff (Elsevier, Amsterdam, 1986). [17] C.B. Duke, A. Paton and A. Kahn, J. Vacuum Sci. Technol. A2 (1984) 515. [18] V.E. de Carvalho, M.W. Cook, P.G. Cowell, O.S. Heavens, M. Prutton and S.P. Tear, Vacuum 34 (1984) 893. [19] V.E. de Carvalho, M. Prutton and S.P. Tear, Surface Sci. 184 (1987) 198. [20] D.J. Titterington and C.G. Kinniburgh, Computer Phys. Commun. 20 (1980) 237. [21] E. Clementi and C. Roetti, At. Data Nucl. Data Tables 14 (1974) 177. [22] J.B. Pendry, J. Phys. C4 (1971) 3095. [23] C.B. Duke, A~ Paton, W.K. Ford, A. Kahn and G. Scott, Phys. Rev. B24 (1981) 3310. [24] C. Mailhiot, C.B. Duke and D.J. Chadi, Surface Sci. 149 (1985) 366. [25] J.B. Pendry, J. Phys. C13 (1980) 937. [26] J. Philip and J. Rundgren, in: Determination of Surface Structures by LEED, Eds. P.M. Marcus and F. Jona (Plenum, New York, 1980) p. 409. [27] A.C. Sobrero and W.H. Weinberg, ref. [26], p. 437. [28] F. Jona, P.M. Marcus and H.D. Shih, ref. [26], p. 247. [29] C.B. Duke, A. Paton, A. Kahn and D.-W. Tu, J. Vacuum Sci. Technol. B2 (1984) 366.