Faraday Discuss., 1993,96, 151-159 Dynamics of Dissociative Chernisorption Cl,/Si( 111)-2 x 1 Michael C. Payne and Ivan Stich Cavendish Laboratory, Madingley Road, Cambridge, UK CB3 OHE Alessandro De Vita and Michael J. Gillan Department of Physics, University of Keele, Keele, Staflordshire, UK ST5 5BG Lyndon J. Clarke Edinburgh Parallel Computer Centre, University of Edinburgh, MayJieEd Road, Edinburgh, UK EH9 3JZ We present a new theoretical methodology for investigating the process of dissociative chemisorption which consists of performing fully ab initio dynamical simulations of the dissociation process. The advantages of the new methodology are: the inclusion of the degrees of freedom of both the incident molecule and the atoms in the surface; an extremely accurate repre- sentation of the interatomic potentials over the relevant regions of configu- ration space and the absence of any assumptions about the nature of the dissociation process.The limitations of the technique are : the restriction of ground-state electronic configurations ; the use of classical dynamics to evolve the ionic trajectories and inadequate statistics. It is shown that this final limitation is already being overcome by the continuing increase in the power of modern computers. The complexity of dissociative chemisorption' -2 and the short timescales sometimes involved in this process make it extremely difficult to investigate the full details of the process using only experimental methods. Theoretical techniques provide an extremely useful complement to experimental techniques for a number of reasons.First, theoretical methods are able to deal with very short timescales although they can themselves run into significant difficulties when a range of timescales are involved. A further advantage is that theoretical systems are extremely well characterised whereas experimental systems inevitably contain a range of different atomic environments. Of course, if the interesting science occurs only at a very small number of exceptional sites the theoretical method is likely to miss the important reaction route until these sites have been identi- fied experimentally. Once this has been done these sites can be incorporated in the theoretical system. In the subject area of this particular meeting, theoretical work has played an important role in both rationalising and motivating experiments.Examples are the calculation of trajectories on potential-energy surfaces (PES)3 and dynamical simulations of molecule-surface collisions performed using model interatomic poten- tial~.~As with most theoretical and experimental techniques there are systems for which these techniques work extremely well and others for which they are not so well suited. In the present paper we describe an alternative theoretical methodology that has its own limitations but which complements these existing techniques. Before moving on to describe this methodology we shall first briefly describe the strengths and weaknesses of existing theoretical techniques.151 Dissociative Chemisorption Cl,/Si( 11 1)-2 x 1 In the PES approach a molecule-surface collision is investigated by evolving a quantum mechanical wavefunction Y({RJ),where R, are the positions of the atoms in the system, according to a time-dependent Schrodinger equation in which the potential energy operator VP({R,})is the total energy of the system with atoms in positions {R,) and the electronic configuration in its Pth eigenstate. It is the term Vp({R,})in the Schrodinger equation which is the PES. One restriction of this approach is immediately apparent, the dynamical evolution of the ionic configuration occurs with the electronic system in a particular eigenstate. Therefore, no electronic transitions and hence no loss of energy via electronic transitions are included in this approach.Energy loss to phonons is included in this technique since the PES does describe these excitations in the surface. However, as the problem stands the PES has a dimensionality of 3N, where N is the total number of atoms in the system including the molecule and the surface. It is well known that when performing an integral in a large number of dimensions, once the dimensionality of a function becomes large it is mathematically intractable to perform the integration numerically by determining the function uniformly in the multi-dimensional space. The same difficulty is present in the PES approach. For instance, even if the PES could be constructed by a knowledge of the total energy for just three points in each spatial dimension for each atom in the system, a hopelessly unrealistic underestimate, then this would still require the calculation of 2 x total energies for a 10 atom system! It is clear that the PES approach is only tractable with severe trunca- tion of the number of degrees of freedom in the system and this implies that only infor- mation about a few of the atoms in the system can be included when constructing the PES.Hence, the PES-based approach is usually applied to diatomic molecules and even then the full six-dimensional PES is often simplified by considering lower dimensional slices through the surface. k consequence of this reduction in dimensionality is that the positions of the surface atoms do not change during the molecule surface collision and excitation of phonons is not included in the formalism.This approximation becomes exact in the limit where the ratio of the mass of the atoms in the molecule to the mass of the atoms in the surface tends to zero. Hence, the truncated PES approach is valid in the case of molecules formed from light atoms hitting surfaces composed of heavy atoms and it is not surprising that the vast majority of the work carried out with the PES approach is for surface collisions involving hydrogen molecules. In the example of performing integrals in large numbers of dimensions the solution to the problem encountered with the uniform-grid sampling method is well known. Rather than attempt to sample the function regularly throughout space, the function is sampled at random points in the function space and the average of the function at these sampling points is used to represent the integral. This is the basis of the Monte Carlo sampling technique which provides an extremely efficient, mathematically well behaved technique for evaluating integrals of multi-dimensional functions.It should be stressed that additional information about the nature of the function can be incorporated into the sampling method in order to increase the efficiency of the Monte Carlo technique, this is generally referred to as importance sampling. Furthermore, information gathered during the sampling can be used as the basis of importance sampling in order to increase the accuracy achieved in later sampling steps. However, in the absence of any such information the function is sampled at random points in the function space.In molecule-surface collisions where the surface cannot be regarded as static the solution to the intractably high dimensionality of the PES is also based on a random sampling procedure as used in the Monte Carlo method described above. Instead of calculating the PES at a regular grid of points the PES is essentially sampled by per- forming a set of trajectory simulations from random initial configurations. In each tra- jectory the system evolves along a line in the muti-dimensional space of the full PES, the crucial point being that only accessible points in the multi-dimensional space are sampled using this technique and inaccessible points are never sampled.In principle, M. C. Payne et al. trajectories should be generated using this approach until the statistical probabilities of the various outcomes of the molecule-surface collision are not changed by simulations of further trajectories. In practice, an importance sampling technique will almost certain- ly be adopted and the procedure will be to attempt to determine which ranges of initial configurations give rise to the same outcome. This can be done by searching for the critical points in the space of initial configurations, where the outcome of the molecule- surface collision changes. A significant reduction in the amount of information required by this approach, compared to a calculation of the full PES, only occurs for classical systems since information at any time-step during the trajectory is then only required for a single position of each ion core, rather than the range of positions required if the ion core is described by a wavepacket.Therefore, this approach can only be adopted for classical ion dynamics and hence it is only meaningful if the quantised energies of the vibrational modes in the system are much smaller than the energy associated with the dissociation process. This approximation becomes more valid as the masses of the atoms in the system increase. This condition is not exactly the inverse of the previous condition for the validity of the PES-based approach since there it was the ratio of the masses of the atoms in the molecule to the masses of the atoms in the surface that provided the criterion for validity of the method. However, there is a degree of complementarity in these two sets of conditions which means that the use of one or other approach allows many more systems to be studied than is possible if just one of these approaches is available. The arguments presented above show that for any molecule-surface collision in which the surface is not static the only tractable approach to studying the collision process is to calculate a series of classical dynamical trajectories from random starting configurations initially until sufficient information has been amassed to concentrate on the most important ranges of initial configurations.Therefore, the calculation of a set of trajectories for various initial conditions for the molecule must be viewed as the correct mathematical method of dealing with the complexity of the problem rather than just a means of producing pictures of molecule-surface collisions. The previous discussion does not depend on the nature of the interactions between the atoms in the system. In the PES approach, these interactions are hidden in the term VP((R,))(i.e. the PES) and in the dynamical approach these interactions determine the forces between the atoms in the system which, in turn, govern the evolution of the atomic coordinates. What has not been stated above is the difficulty associated with determining these interactions. There are a number of techniques for calculating inter- atomic interactions ranging from full ab initio calculations to simple ball and spring models.While it is true that any correctly parametrised model of the interatomic inter- actions is correct in some region of phase space, the difficulty is that the extent of this region of validity is not known. In a process such as dissociative chemisorption signifi- cant changes in bond-lengths must occur and the environment of the atoms in the mol- ecule changes completely during the dissociation process. In order to describe the dissociation process correctly, therefore, the interatomic potential must be correct over a very wide region of configuration space. This places extreme demands on the inter- atomic potential to be used in the simulations.Without entering the debate about the validity of various classes of interatomic potential, there is one crucial point to be made. If it is decided to run trajectory simulations with a particular model potential then it may be necessary to verify the correctness of this model potential at each of the critical points during the reaction process if the configurations at these critical points are not clearly within the range of validity of the model potential. In order to perform this verification a calculation must be performed with a more accurate and reliable method. If the model potential is found to be incorrect at any of these critical points the entire trajectory must be rejected. Therefore, even if a model potential is used to increase the speed of the simulation, it is possible that full verification of the findings could incur as Dissociative Chemisorption CI,/Si( 1 11)-2 x 1 much cost as performing the entire simulation using a more accurate and reliable poten- tial.Of course, this final verification step is rarely performed. If this is the case, the trajectory simulations using the simple model potential are very fast computationally but there may be doubts about the validity of the results. These problems do depend on the nature of the system and there are many systems for which simple model potentials work very well. Particularly noteworthy are the various ‘glue’ model potentials for metallic systems which are remarkably successful at describing the bonding over an enormous range of configuration space.’ In contrast, there is at present no simple model potential for the tetrahedral semiconductors that works well over a large region of con- figuration space and it is generally necessary to use ah initio methods to determine quantitative values of many physical properties in these systems.The previous argument suggests that there is a good reason for using the most reli- able technique for calculating interatomic potentials for trajectory simulations. The argument against this is not so much that of the increased computational cost of, say, ah initio methods, but more that such methods cannot be applied to systems that are large enough to model the dissociation process.For instance the configuration interaction method, which can be exact, may only be applied to a few-electron system before absol- ute accuracy is lost. Less costly computational techniques such as the Hartree-Fock method can be applied to much larger numbers of electrons and provide very good accuracy for closed-shell systems but are relatively poor for other systems. In the longer term, the favourable scaling of quantum Monte Carlo methods may allow essentially exact calculations to be performed for large systems. It should also be mentioned that physical properties calculated with quantum chemistry methods which use a cluster to represent a surface do not converge rapidly with cluster size. The convergence can be significantly improved by careful termination of the cluster, in the case of silicon this is very effectively achieved by saturating the dangling bonds of the surface of the cluster with hydrogen atoms.Computational methods that have developed from solid-state physics do not, at first sight, appear to have any relevance to the problem of dissociative chemisorption since these techniques can generally only be applied to systems which have full three-dimensional periodicity. However, by performing calculations on a periodic sequence of slabs with vacuum gaps in between it is possible to create a set of surfaces in a system which has full three-dimensional periodicity and if a molecule is placed in the vacuum gap moving towards the surface all the features required to simu- late dissociative chemisorption are present. Of course, creating such a geometry within the constraint of the three-dimensional periodicity entails an enormous computational cost since the unit cell must contain at least a hundred atoms and a large vacuum region.Such a system size is presently beyond most ah initio methods. However, one computational technique has developed so rapidly during the last ten years that it can be applied to systems of this size. This method is the total energy pseudopotential tech- nique.6 The features of this technique that make it computationally efficient are: (i) the use of pseudopotentials to represent the ions so that only the wavefunctions of the valence electrons have to be computed; (ii) the use of density functional theory to map the interacting system of electrons onto a system of non-interacting electrons moving in an effective potential, V,, ,the exchange-correlation potential.This mapping replaces the direct product wavefunction of the true electronic system by a simple product wavefunc- tion, hence the dimensionality of the wavefunction space is drastically reduced and, more importantly, scales linearly with the number of electrons in the system. In principle this mapping is exact and density functional calculations should give the exact ground- state energy and electronic density. In practice, the form of the exchange-correlation potential is not known exactly and so approximations have to be adopted. The approx- iations typically give physical quantities to an accuracy of a few per cent.One restriction incurred by the use of density functional theory is that the interatomic interactions can only be computed for the ground-state electronic configuration, hence only one of the M. C. Payne et al. many PES can be computed; (iii) the use of a plane wave basis set to expand the elec- tronic wavefunctions. A plane wave basis set is possible because of the use of pseudo- potentials and the three-dimensional periodicity of the system. Although a large number of plane waves are needed for each atom in the system the numerical operations on these basis functions are extremely fast and computationally efficient. Hence, the overall computational cost of the total energy pseudopotential method scales very favourably with system size.Furthermore, it is very easy to calculate the forces on the ions with a plane wave basis set since there are no Pulay forces; (iv) the ability to implement total energy pseudopotential calculations scalably on parallel computers. This is important since there is then, in principle, no computational limit to the size of system that can be studied using this method. In practice, since bigger systems require bigger computers there is effectively a cost limit to the size of accessible system but even present gener- ation parallel machines are powerful enough to perform calculations on systems con- taining several thousand atoms which places an enormous number of scientific problems, including many surface science problems, within the scope of the method.This extended introduction has demonstrated the need for trajectory-based classical ionic dynamics simulations for studying dissociative chemisorption in a system in which there is a significant surface response and the need to perform high quality calculations of interatomic interactions in order to have faith in the results of these simulations. For these reasons we have chosen to investigate the ability of the total energy pseudo- potential method of modelling such a complex process by performing a number of fully ab initio dynamical simulations of chlorine molecules striking the 2 x 1 n-bonded (1 11) surface of silicon.’ Simulations The previous section presented the case for ab initio simulations of dissociative chemi- sorption but the enormous computational cost associated with such a study must be considered.Our calculations were performed on a 64-node Meiko Computing Surface and a 64-node Intel Hypercube. Each trajectory required of the order of 10 days com- puting at a rate of 1 GFlop. These may appear to be enormously large numbers and might be expected to be associated with an astromonically high cost. However, if an attempt is made to quantify the cost per trajectory this can be estimated (determined by the total cost of the Meiko machine divided by the proportion of its three year oper- ational life expended on each trajectory simulation) to be of the order of ElOOOO. We shall consider this further in the following section. Given the limited computational time available to us we were able to carry out just five trajectory simulations.Thus it was imperative to ensure a high probability of dissociation in each simulation. Therefore, our choise of system and incident energy of 1 eV for the chlorine molecule was motivated on the basis of experimental results that showed a very high probability of dissociation at this energy in related sy~tems,~~~ Furthermore, at this incident energy the dissociation process must be chemically driven and so a correct description of the electronic struc- ture during the dissociation process is expected to be important. We also chose to study a semiconductor surface since our formalism can only describe the electronic ground- state configuration and low-energy electron-hole pairs cannot be generated in this system. We note, however, that no method is presently capable of including electron excitation during the dissociation process.Such a limited number of trajectories makes a mockery of a statistical approach to sampling the multi-dimensional PES. However, the types of questions we can address are: (i) Is the dissociation instantaneous, occurring as the molecule strikes the surface, or does the molecule first get trapped in a precursor state? (ii) Do extremes of initial orien- tation affect the nature of the dissociation process? (iii) What is the mechanism of energy dissipation in the dissociation process? Given the restricted number of trajectories we do Dissociative Chemisorption Cl,/Si( 1 11)-2 x 1 not include rotational or vibrational motion in the incident molecule.The initial orien- tations of the chlorine molecule for our five trajectories are shown in Table 1. Snapshots from the five trajectories, giving simple ball and stick representations of the atomic positions, are shown in Fig. 1. For each trajectory, Fig. 1 shows the initial configuration, an intermediate configuration and the final configuration generated in the simulation. Since the purpose of the present paper is to motivate this new theoretical methodology only a few comments about the results of the simulations will be made here. The chlo- rine molecule dissociates in each of our simulations; this is not surprising given that we have chosen an incident energy at which it is known experimentally that almost 100% dissociation occurs. In trajectories 1, 2 and 3 the dissociation of the chlorine molecule is instantaneous but in trajectories 4 and 5 the molecule first binds to the surface, under- going significant extension of the molecular bond, and dissociation only occurs when the molecular axis rotates towards the surface.In trajectory 1 the chlorine atom that moves across the surface is captured only as a result of the motion of the silicon atom to which it is eventually bound. If the surface was assumed to be static this chlorine atom would have escaped back into the vacuum. In all the simulations the effect of rehybridisation can be clearly seen; the silicon atoms bound to the chlorine atoms change from their almost planar configuration in the original 2 x 1 reconstructured surface to a tetra- hedral configuration as a result of the change from sp2 to sp3 hybridisation.In trajec- tory 5 it can be seen that the impact of the chlorine molecule with the silicon surface induces a change in the bonding topology. As mentioned previously, our simulations are performed using three-dimensional supercells. This means that the system actually studied is a periodic array of slabs and a three-dimensional periodic array of chlorine molecules are incident on the surfaces of these slabs. We also impose inversion symmetry on the system so that chlorine mol- ecules are simultaneously incident from either side of the slab.Our supercell dimensions are 11.5 A x 13.3 in the plane of the surface and 24.5 A normal to the surface. The slab consists of four double layers giving a total of 96 silicon atoms in the unit cell. The vacuum gap between successive slabs is 14.3 A. Test simulations performed with smaller systems showed effects which could be attributed to the periodic boundary conditions but no such effects were observed in the large systems used for our simulations. The simulations were performed at constant total energy so that the energy released by the dissociation of the molecule is converted to thermal energy in the slab. This leads to an increase in the temperature of the system during the simulation of several hundred K. The trajectories were followed for between 200 and 400 fs and were terminated when the outcome of the reaction was known, generally when both chlorine atoms were bound to the surface.Finally, we note that by time-reversing our trajectories we describe the recombination of chlorine atoms on the silicon surface to chlorine molecules in the vacuum. Table 1 Initial position of chlorine molecule relative to the silicon surface trajectory 1 2 3 4 5 initial position of C:P:I V:P:I V:P:II H:N C:N chlorine molecule The first letter denotes whether the chlorine molecule is above a n-bonded chain (C), a valley between the n-bonded chains (V), or a hole formed by the six- membered ring that extends from the chain to the valley (H). The second letter denotes whether the chlorine molecule is parallel (P) or normal (N) to the surface.For the configurations in which the chlorine molecule is parallel to the surface the final symbol denotes whether the molecule is parallel (11) or perpendicular (I)to the n-bonded chains M. C. Payne et al. Fig. 1 Ball and stick representations of the starting, intermediate and final atomic configurations in our five simulations. The trajectories are shown in the same order as in Table 1. Dissociative Chemisorption Cl,/Si( 11 1)-2 x 1 One reaction to our work is that all the results of the simulations are totally predict- able. Given the impossibility of predicting even the small amount of information pre- sented in Fig. 1 we find this claim unbelievable. The claim is easily stated after the simulations have been performed and it is indeed reassuring that the results can be rationalised according to sound chemical and physical principles but that is not the same thing as predicting the results before the simulations were carried out.The simula- tions reveal phenomena which depend critically on both the atomic and electronic response of the surface response which would not have been discovered using any alter- native methodology. Of course, the system was chosen purely because existing theoreti- cal techniques could not be expected to provide an accurate description of the dissociation process. The authors have never intended the ab initio dynamical approach to displace existing techniques, all of which have particular strengths and weaknesses and work well for some systems and poorly for others.This is just as true of our method as of any other method. Indeed, we have ensured that our computational technology has been made available so that realistic PESs can be computed for systems for which the PES approach is the correct methodology." However, there are systems for which there is at present no other viable theoretical methodology for studying dissociative chemi- sorption and for these systems there should be no argument about the usefulness of our approach. If our results were so predictable we would like to see a paper published before our calculations were carried out that predicts the results of the simulations. Furthermore, if this system is so trivial and predictable then it must be equally simple to deal with the case of 02/Si(lll)-2 x 1 so we would like to see a complete description of the process of dissociative chemisorption for 02/Si( 11 1)-2 x 1.The Future The enormous advances in the total energy pseudopotential method have arisen from a combination of algorithmic improvements, initiated by Car and Parrinello,' ' and an increase in the power of computers, most notably through the development of parallel computers. The advance in computer power has not occurred uniformly across the spec- trum of machines. The speed of expensive supercomputer processors has not increased significantly in the last decade but during this time the speed of cheaper processors has become almost equal to that of supercomputer processors.These cheaper processors thus have an enormous price to performance advantage over supercomputer processors. Parallel computer manufacturers are able to exploit this performance advantage and assemble large number of processors at modest cost into a machine that significantly outperforms conventional supercomputers. The same price to performance advantage means that it is now possible to buy a workstation for a modest sum of money that is as powerful as a supercomputer. If we were to perform more trajectory simulations using the most recently released workstations the cost per trajectory (calculated on the same basis as the previous figure) would be of the order of E1000. This reflects the fact that in just three years the cost of computing has decreased by a factor of 10.This rate of decrease has been roughly the same for the last decade. Although there must be some slowing of this trend eventually, if it continued at the present rate then the cost per trajectory would be &lo0in three years time, El0 in six years time and El in nine years time. Even if the rate of decrease in the cost of computing slows somewhat the cost of 510 per trajectory will be reached relatively soon. It should be clear that the most serious limitation of our simulations, that of non-existent statistics, is not a long-term problem and that the cost effectiveness of ab initio theoretical simulations compared to experiment will not only increase but increase very rapidly.Our calculations verify the theoretical methodology and prove that the process of dissociative chemisorption is accessible to ab initio simulation. This is important in itself, but the real excitement is the certainty that so much more can be done in the future. Furthermore, there will be more M. C. Payne et al. and more scientific problems, including many involving surface dynamics, where theo- retical calculations will eventually be more cost effective than experiment. The simulations were performed as part of the ‘Grand Challenge’ collaborative project, coordinated by Prof. D. J. Tildesley, on the Meiko Computing Surface at Edinburgh University. Some calculations were also performed on the Intel iPSC/860 Hypercube at the Daresbury Laboratory.We acknowledge financial support from the SERC under grants GR/G32779 and GR/E92582. We thank Dr. K. C. Bowler, Dr. S. P. Booth and Dr. P. J. Durham for their support. References 1 See for instance, Interactions of Atoms and Molecules with Solid Surfaces, ed. V. Bortolani, N. H. March and M. P. Tosi, Plenum Press, New York, 1990. 2 For a recent review see A. W. Kleyn, J. Phys.: Condens. Matter, 1992,4, 8375. 3 S. Holloway, J. Phys.: Condens. Matter, 1991, 3, S43. 4 F. H. Stillinger and T. A. Weber, Phys. Rev. Lett., 1989,62,2144. 5 I. J. Robertson, M. C. Payne and V. Heine, Europhys. Lett., 1991, 15, 301. 6 M. C. Payne, M. P. Teter, D. C. Allan, T. A. Arias and J. D. Joannopoulos, Rev. Mod. Phys., 1992, 64, 1045. 7 K. C. Pandey, Phys. Rev. Lett., 1981,47, 1913. 8 M. L. Yu and B. N. Eldridge, Phys. Rev. Lett., 1987,58, 1691. 9 D. A. Hansen and J. B. Hudson, Surf Sci. Lett., 1991,254,222. 10 B. Hammer, K. W. Jacobsen and J. Norskov, Phys. Rev. Lett., 1992,69, 1971. 11 R. Car and M. Parrinello, Phys. Rev. Lett., 1985,55, 2471. Paper 3/03008C; Received 24th May, 1993
展开▼
