Faraday Discuss., 1997, 106, 135»154 Ab initio simulation of molecular processes on oxide surfaces Philip J. D. Lindan,a* Joseph Muscat,a Simon Bates,b Nicholas M. Harrisona and Mike Gillanb a DCI, CCL RC Daresbury L aboratory, W arrington, UK WA4 4AD b Physics Department, Keele University, StaÜordshire, UK ST 5 5BG Ab initio calculations based on both density functional theory (DFT) and Hartree»Fock (HF) methods are used to investigate the energetics and equilibrium structure of the stoichiometric and reduced (110) surface, the TiO2 adsorption of potassium on the (100) surface and of water on the (110) surface.It is shown that DFT and HF predictions of the relaxed ionic positions at the stoichiometric surface agree well with each other and fairly well with recent X-ray diÜraction measurements.The inclusion of spin polarisation is shown to have a major eÜect on the energetics of the reduced surface formed by removal of bridging oxygens. The gap states observed to be induced by reduction are not reproduced unless spin polarisation is included. Static and dynamic DFT calculations on adsorbed water on TiO2 (110) con–rm that dissociation of leads to stabilisation at low cover- H2O ages, but suggest a more complex picture at monolayer coverage, in which there is a rather delicate balance between molecular and dissociated geometries.I Introduction Until recently, the understanding of oxide surfaces lagged far behind that of metals and semiconductors. This was partly because of experimental difficulties, and partly because the complexity of many oxide materials made them difficult to study with accurate theoretical methods.The situation is now changing rapidly, and the last few years have seen important new experimental measurements1h3 as well as accurate quantum-mechanical calculations on a variety of oxide surfaces and on adsorbed molecules.4 Titanium dioxide is of particular interest, both because of its important and varied applications, and because of the challenges it poses to ab initio theory. The theoretical problems arise from the treatment of exchange and correlation in partially occupied d states when the material is reduced.These problems make it interesting to compare diÜerent ab initio methods. We have therefore used both DFT and HF approaches using plane-wave and local basis sets in order to explore the sensitivity of the computed answers to the approximations made. Our calculations address three issues : the equilibrium structure of the perfect surface, the atomic and electronic structure of the reduced surface obtained both by removal of oxygen and by addition of alkali metal and the energetics and structure of water adsorbed on the surface.Early ab initio calculations on the energetics and relaxed structure of oxide surfaces were made on MgO (001) using HF theory,5 and HF calculations have since been reported on several other oxide surfaces, including The HF approach has TiO2 .6h10 also been used in ab initio studies of molecular adsorption, work on adsorption at H2O MgO (001)11 being particularly relevant to the results to be reported here.DFT methods began to be applied to oxide surfaces somewhat later, but progress has been rapid, and a number of detailed studies have been reported on perfect and defective oxide surfaces, as well as on the adsorption of simple molecules.12h23 135136 Simulation of molecular processes on oxide surfaces The equilibrium positions of the ions at oxide surfaces can diÜer greatly from their perfect-lattice positions, and the comparison of predicted and measured positions is an important way of testing the reliability of the theoretical methods. Recent progress in the application of glancing angle X-ray diÜraction to oxide surfaces provides detailed surface structural information for the –rst time. We shall compare DFT and HF predictions of the relaxed positions at the perfect (110) surface with the very recent TiO2 diÜraction measurements of Charlton et al.3 The comparisons demonstrate both the good accord of the two theoretical approaches and their semi-quantitative agreement with experiment.Departure from stoichiometry is a characteristic of transition-metal oxides. Reduction by removal of oxygen is particularly important, since experiments are often performed under ultra-high vacuum, which tends to induce oxygen loss.In addition, many experiments require signi–cant surface conductivity, which is often produced by deliberate reduction. Reduction is also important in practice, since is often sub- TiO2 stoichiometric in real applications or under processing conditions. We shall present DFT and HF calculations on the relaxed structure of the reduced (110) surface TiO2 and on its electronic structure.A major issue here is that of spin polarisation, since reduction from Ti4` to Ti3` necessarily leads to unpaired electrons, and we shall show that the inclusion of spin polarisation in the calculations is essential. Reduction can be produced not only by oxygen removal but also by the adsorption of reducing species such as alkali metals, and we shall present HF calculations on the eÜects of K adsorption on the (110) surface.TiO2 Finally, we shall address the problem of water adsorption. Recent static20 and dynamical calculations22 indicate that dissociates rather readily on (110). We H2O TiO2 present new calculations which con–rm this, but which also suggest that the balance between molecular and dissociative adsorption is quite delicate.II Techniques The DFT-pseudopotential and HF techniques for calculating the energetics of solids and their surfaces have been reviewed extensively,24,25 and here we give only a brief summary. Both techniques work with single-electron orbitals which satisfy a ti(r) Schroé dinger-like equation: Hti\A[ +2 2m Z2]Ven]VH]V effBti\viti (1) where the potential operator is separated into the electron»nucleus interaction the Ven , Hartree potential and an operator representing the eÜects of exchange and VH , Veff correlation. However, the meanings of the orbitals and the operator are diÜerent ti Veff in the two approaches.In DFT, the are not real single-electron wavefunctions, but are merely auxiliary ti(r) quantities used to represent the electron density n(r).The exchange-correlation operator has the form of a local potential which is de–ned to be the functional derivative Veff(r) where is the exchange-correlation energy. The major issue in DFT is the dExc/dn(r), Exc approximation used for Exc . Until recently, most DFT work on solids was based on the local-density approximation (LDA):26 Exc\Pdr n(r)vxc[n(r)] (2) where is the exchange-correlation energy per electron in a uniform electron gas of vxc(n) density n.This works well for many bulk solids, but is not accurate enough for molecular dissociation energies, and most of the work to be described later is based on gener-P. J. D. L indan et al. 137 alised gradient approximations (GGAs),27,28 in which the dependence of on Exc gradients of n(r) is included. A widely used implementation of DFT combines a plane-wave basis set with the pseudopotential method.29,30 The pseudopotential method assumes that the core electrons are in the same states as in the free atoms, and the pseudopotential represents the interaction between valence electrons and the ionic cores. The pseudopotential is generated via DFT calculations on free atoms.Errors can arise from the choice of pseudopotential, but these errors are usually under satisfactory control. In the pseudopotential approach, it is common practice to use supercell geometry (i.e. periodic boundary conditions) and plane-wave basis sets, so that the occupied orbitals are expanded ti(r) as : ti(r)\; G ciG exp(iG … r) (3) where the sum goes over reciprocal lattice vectors G of the supercell lattice.This is an in–nite sum in principle, but in practice it is restricted to those G for which +2G2/2m\ where is a chosen cut-oÜ energy. The size of governs the completeness of Ecut , Ecut Ecut the basis set. The practical determination of the ground state is performed by minimising the total energy of the system with respect to the plane-wave coefficients ciG .For some of the work on it is essential to include electronic spin. The gener- TiO2 , alisation of DFT to include spin polarisation is formally straightforward, and GGAs in which depends on the densities and of up-and down-spin electrons are well Exc nè(r) né(r) established.In practice, we perform the calculations by –xing the separate numbers of spin-up and spin-down electrons in advance. Comparison of ground-state solutions with diÜerent spin excesses is then made to determine the lowest-energy solution. In contrast to DFT, the HF approximation attempts to describe the many-electron wavefunction . . . , The orbitals appearing in eqn. (1) are the factors in the W(r1, rN).ti antisymmetrized product representing W: W(r1, . . . , rN)\AA< i/1 N ti(ri)B (4) The operator now represents the non-local exchange interaction between single- Veff electron orbitals : Veffti(r1)\Pdr2 V (r1, r2)ti(r2) \[;j Pdr2 e2 o r1[r2 o tj p (r2)tj(r1)ti(r2) (5) By de–nition, electron correlation is neglected in this approach, though approximations are available for reintroducing correlation.31h33 It is common practice to represent the in a local basis set of Gaussian orbitals, in which all the required integrals ti are analytic.Spin-polarised calculations are readily performed in the HF approximation. Because the energy expression is based on analytic integrals, the numerical stability is very high, and very small magnetic coupling energies can be studied.DiÜerent spin states can be explored by initialising the calculation with a superposition of spinpolarised ions. Our DFT-pseudopotential calculations have been performed using two separate codes: the CASTEP code24 (or its parallel version CETEP34) and the VASP code.35 The two codes employ essentially the same overall strategy, namely a ground-state search by global minimisation of the total energy with respect to the plane-wave coefficients.The diÜerences are that VASP uses more efficient search techniques, and is also able to use Vanderbilt-type ultrasoft pseudopotentials,36 which means that a given accuracy can be138 Simulation of molecular processes on oxide surfaces Fig. 1 (110) surface of Light and dark spheres indicate titanium and oxygen ions, respec- TiO2 .tively. (a) Perspective view showing the slab geometry used. The simulation cell is extended for display purposes. The –ve-fold and six-fold coordinated titanium sites and the bridging oxygen site are labelled 5f, 6f and BO, respectively. (b) The 1]1 surface unit cell used in the simulations. (c) Side view of the relaxed stoichiometric geometry. achieved with much smaller plane-wave basis sets.The HF calculations described later were carried out using the widely used CRYSTAL package,25,37 which is also capable of performing DFT calculations. III TiO (110) surface Several research groups have recently used ab initio methods to investigate the relaxed structure of the (110) surface. Surface X-ray diÜraction measurements of the TiO2P.J. D. L indan et al. 139 relaxed ionic positions at the surface have given crucial evidence that the theoretical predictions are essentially correct. We compare here the experimental values of ionic displacements with published DFT predictions,17 and with our own previously unpublished HF and DFT-pseudopotential predictions. The structures of the perfect crystal and of the perfect stoichiometric surface are shown in Fig. 1. A prominent feature of the surface is the rows of ìbridgingœ oxygens which stand proud of the surface plane of Ti and O ions. The latter plane contains two kinds of Ti ions : six-fold coordinated Ti ions lying beneath bridging oxygens, and exposed –ve-fold Ti ions ; the O ions in this plane are referred to as in-plane O.Experiment and previous calculations show substantial displacements of these surface ions away from their perfect lattice positions, and they also show that six-fold Ti and inplane O move out of the surface, while –ve-fold Ti and bridging O ions move in. Our new DFT-pseudopotential calculations were performed with Vanderbilt ultrasoft pseudopotentials using the VASP code.35 Exchange and correlation were included using the Perdew-91 GGA functional.27 We employed the usual supercell method, and the surface was treated using periodically repeated slab geometry (see Fig. 1). In this scheme, it is absolutely crucial to ensure that the slab thickness, L , and the width, L @, of the vacuum layer between neighbouring slabs are large enough so that interactions between surfaces are negligible. We have systematically studied the convergence of the equilibrium surface structure as L and L @ are increased.We –nd that convergence with respect to L @ is very rapid, and a vacuum width of 4 reduces errors in relaxed posi- ” tions to ca. 10~2 Convergence with respect to L is more complex, and we –nd oscil- ”. latory eÜects due to interactions through the body of the slab.However, when L is increased beyond ca. 20 the errors in relaxed positions again fall below ca. 10~2 ”, ”. A detailed theory»experiment comparison of the ionic displacements is given in Table 1. We note the following points : –rst, independent DFT-pseudopotential calculations are in quite satisfactory agreement with each other ; second, DFT and HF calculations are also in reasonable accord; third, the theoretical and experimental displacements are in semi-quantitative agreement, except for the bridging oxygens.The reason for the substantial discrepancy for bridging oxygens is unclear, but it is evidently signi–cant, since all the calculations show roughly the same discrepancy. Table 1 Ionic displacements due to relaxation of the (110) 1]1 surface (a) (b) (c) (d) label [110] [16 10] [110] [16 10] [110] [16 10] [110] [1 610] 1 0.13 0.00 0.23 0.00 0.13 0.00 0.12 ^0.05 0.00 » 2 [0.17 0.00 [0.11 0.00 [0.13 0.00 [0.16 ^0.05 0.00 » 3 [0.06 0.00 [0.02 0.00 [0.11 0.00 [0.27 ^0.08 0.00 » 4 0.13 [0.04 0.18 [0.05 0.09 [0.07 0.05 ^0.05 [0.16 ^0.08 5 0.13 0.04 0.18 0.05 0.09 0.07 0.05 ^0.05 0.16 ^0.08 6 [0.07 0.00 0.03 0.00 [0.04 0.00 0.05 ^0.08 0.00 » 7 0.06 0.00 0.12 0.00 » » 0.07 ^0.04 0.00 » 8 [0.08 0.00 [0.06 0.00 » » [0.09 ^0.04 0.00 » 9 0.02 0.00 0.03 0.00 [0.01 0.00 0.00 ^0.08 0.00 » 10 [0.03 [0.05 0.00 [0.02 » » 0.02 ^0.06 [0.07 ^0.06 11 [0.03 0.05 0.00 0.02 » » 0.02 ^0.06 0.07 ^0.06 12 [0.01 0.00 0.03 0.00 » » [0.09 ^0.08 0.00 » 13 » » 0.00 0.00 » » [0.12 ^0.07 0.00 » Labels refer to Fig. 1(c). The displacements are in and are from the bulk terminated positions. ”, The results are from: (a) DFT calculations by Ramamoorthy et al. ;17 (b) the present DFT calculations ; (c) the present HF calculations and (d) surface X-ray diÜraction experiments.3140 Simulation of molecular processes on oxide surfaces Overall, the present understanding of the equilibrium structure of the stoichiometric surface seems to us fairly satisfactory, and we think it provides an adequate basis for the work on molecular adsorption and surface reduction reported below.IV Reduced TiO2 We now turn to the reduced (110) surface. As explained in the Introduction, TiO2 reduction can be achieved either by removal of oxygen or by adsorption of reducing species such as alkali metals.We begin by considering oxygen removal. Two questions arise when considering a transition-metal oxide which is oxygen de–- cient. The –rst is, what is the defect structure, both physical and electronic, associated with the reduction ? The second, what are the consequences of the reduction for surface structure, reactivity, transport etc. For it turns out that the defects associated with TiO2 reduction are in themselves very complex, and here we shall con–ne ourselves to the –rst of these questions.Speci–cally, we will look at oxygen vacancies in the bulk and on the (110) surface. The examination of computed results from both DFT theory implemented using pseudopotentials and a plane-wave basis set and HF theory using local basis functions allows us to examine the sensitivity of these results to the main theoretical and numerical approximations made.A familiar feature in experimental studies of reduced is the presence of ìTi3`œ ions. The evidence for this description is discussed in ref. 38, and we summarise the key features here. The UV photoemission spectroscopy (UPS) spectra contain a bandgap feature for reduced samples which is interpreted in terms of occupied states formed from Ti(3d) orbitals,38,39 a view supported by the resonant behaviour of UPS across the range of photon energies corresponding to the Ti 3p»3d excitation threshold.The shift of the Ti(2p) core levels seen in X-ray photoelectron spectroscopy (XPS)40 and the lack of surface conductivity also evidence the localised nature of the additional electrons. Regardless of the location of the oxygen vacancies (i.e.in the bulk or at a surface) the energy levels of the excess electrons lie in the upper half of the bulk bandgap, though their position depends on the degree of reduction. The exact nature of these states, and their relation to the oxygen vacancy, have not been established by experiment. There have been previous theoretical studies of the reduced (110) surface of TiO2 .Most of these predict gap states, but there is no accord on the surface defect structure responsible for these states. Wang and Xu41 (tight-binding extended Hué ckel), and Tsukada et al.42 (DV-Xa cluster methods) both –nd the states 0.7 eV below the conduction-band minimum. However, Munnix and Schmeits,43 using a tight-binding model, found that gap states only occurred after the removal of sub-surface oxygen atoms. It is, of course, highly desirable that –rst-principles methods be able to describe the electronic structure of reduced transition-metal oxides, being a particular case.The TiO2 fact that the Ti4` ion in rutile is in the d0 state, and that upon moderate reduction Ti3` ions in the d1 state are produced, has important consequences for the use of DFT theory.The failure of the LDA approximation to produce the correct ground state in many transition-metal oxide systems is now well known.44 For example, is described La2CuO4 as a non-magnetic metal when in reality it is a magnetically ordered wide band-gap insulator. Although for many years this failure was attributed to an incorrect treatment of electron correlation we now know that a qualitatively correct ground state is obtained if the description of the exchange interaction is improved.45h47 This understanding underpins the alternative corrections to the LDA such as the self-interaction correction (SIC-LDA)46,47 and the inclusion of explicit on-site potentials (LDA]U).45 In the HF approximation, the exchange interactions are treated exactly, and this is therefore of great interest in these systems.48h50 One might expect such corrections to beP.J. D. L indan et al. 141 rather small for the d1 state expected in reduced because the orbital is not exces- TiO2 sively contracted (W. Temmerman, personal communication). It is, therefore, interesting to compare the results of LSDA calculations and HF calculations in some detail.As will become clear from our results, a qualitatively correct description of excess-electron gap states in can be obtained from calculations performed within either HF or DFT approximations. IV. A Oxygen vacancy in bulk TiO2 It is natural to approach the question of reduced through calculations on the bulk TiO2 material, but this immediately raises the question of the defect structure of This TiO2~x .is a surprisingly complex issue, since a family of stable Magneç li phases51 is found to occur between and It is beyond the scope of this work to calculate the Ti2O3 TiO2 . most stable structure for a given stoichiometry, and this is not necessary for the present purposes. Here, we investigate the predictions of spin-polarised DFT and HF for the electronic structure of the reduced bulk, the general features of which should not depend on the exact defect structure.One caveat is that the defect concentrations accessible to simulation are much higher than those encountered in experiment, a fact which should be borne in mind when comparisons are made. Using plane-wave pseudopotential DFT we have treated systems consisting of one, two and eight unit cells, each containing a single vacancy.23 In each case we –nd that the lowest energy state is spin-polarised, with the two excess electrons occupying states formed from Ti(3d) orbitals.For all three defect densities studied, these states are localised on the three Ti ions nearest to the oxygen vacancy. The energies of these states lie in the gap between occupied and unoccupied levels : for the largest system they were 1.3» 1.9 eV above the valence band maximum, and occupied similar positions in the other systems.This behaviour is in contrast to spin-paired calculations on reduced TiO2 . These do not give gap states, and instead yield metallic solutions in which the excess electron states, although still formed from Ti (3d) orbitals, remain in the conduction band, and are not localised.The energy gained through spin polarisation (i.e. the energy drop in going from spin-paired to spin-polarised solutions) ranged from 0.2 to 0.5 eV. For the two unit cell system we performed full structural relaxation. This did not change the character of the excess-electron states, but did raise their energies in the gap slightly.Unrestricted HF calculations on the reduced bulk have also been made, using a system of four unit cells formed by doubling in the M001N plane.10 To make explicit allowance for the possibility that the excess electron density resulting from oxygen removal might localise, at least in part, at the vacancy site, as it does at the MgO (001) surface,12,52 we have carried out calculations using basis sets both with and without the complete set of oxygen functions at the vacant (oxygen) site.For the bulk vacancy, i.e. both basis sets lead to insulating, triplet-spin ground states, with that derived Ti8O15 , from the larger basis set (vacancy functions included) lower in energy by 1.4 eV. The Mulliken population analysis of bulk indicates site charges of Ti(]2.8) and TiO2 O([1.4) and thus, when an O vacancy is created, an excess electron density of [1.4 is generated.In this state [0.8 of the excess electron density is localised at the vacancy in a diÜuse s-like orbital, [0.4 at the nearer of the two Ti sites and [0.1 at the other Ti, while the net spin populations at the three sites are 0.8, 1.0 and 0.2, respectively. Both solutions are triplet-spin polarised, with gap states 1»2 eV above the valence band upper edge.States of this kind were only found in spin-polarised (unrestricted) solutions. The reason for the diÜerences between spin-paired (restricted) and spin-polarised (unrestricted) solutions lies in the nature of the states occupied by the excess electrons.The removal of a neutral oxygen ion leaves two electrons which previously occupied O(2p) levels in the valence band. These states are no longer available, and the electrons must go into the conduction band, the bottom of which is formed from Ti(3d) orbitals.142 Simulation of molecular processes on oxide surfaces Without spin polarisation the excess-electron states are forced to be doubly occupied.The spin-polarised solution also has the excess electrons occupying Ti(3d) orbitals, but now these electrons are unpaired, which is the origin of the large energy diÜerence between the two solutions. The occupation of spin-unpaired orbitals rather than the partial occupancy of spin-paired orbitals reduces the on-site exchange energy but increases the kinetic energy.In these narrow d-bands it is clear that the exchange term dominates and, therefore, the ground state has spin-polarised electrons which are localised on the Ti sites.53 The energy gained lowers these states. A further consequence is that the spin-unpaired orbitals contract closer to the ionic core, which is apparent through comparison of the excess-electron charge density distributions from the spinpaired and spin-polarised calculations (not shown here).Whilst yielding results that are qualitatively similar, there are notable diÜerences between the DFT and HF solutions for the reduced bulk. HF predicts that at least some of the excess-electron charge occupies the oxygen vacancy site, whilst in all three DFT calculations, negligibly small charge density was located there.This diÜerence indicates that the nature of the electronic state is sensitive to the approximations made. We are currently examining the cause, which may be due to the diÜerent one electron potentials or to the numerical approximations (basis set, pseudopotential). The variation in the position of the excess-electron levels above the valence band is of less concern: the signi–cant points are that in both cases they are separated from unoccupied levels, and that this separation arises from spin-polarisation.IV. B Reduced (110) surface The surface calculations we report here have used periodically repeating slab geometry, as indicated in Fig. 1 and discussed previously. The reduced 1]1 surface is formed by removing all the bridging oxygen ions, giving a density of surface vacancies of one monolayer, and a total of four excess electrons in the slab.We have performed both DFT23 and HF10 calculations on this surface, and in both studies we investigated solutions with spin-polarisation and full structural relaxation. Here, we make direct comparison of the predictions of the methods.In addition, the plane-wave study includes spin-paired solutions, allowing us to assess the eÜects of spin-polarisation, both at –xed geometry and after full relaxation. Both methodologies yield ground-state solutions that have insulating spin-polarised states above the valence band maximum, and these states localise the excess charge roughly equally on the two surface titanium ions, in bands formed from Ti(3d) orbitals.According to HF, these states lie 2 to 3 eV above the valence band maximum, while DFT places them 0.71 to 1.8 eV above; the lowest unoccupied states are ca. 1.9 eV above the valence-band edge, less than the experimental bandgap (3.1 eV)39 as is usual for DFT eigenvalues. These latter results are illustrated in Fig. 2, which shows the calculated density-of-states for up-and down-spin electrons.The excess-electron states are clearly visible in the up-spin density, and some spin-polarisation of lower levels is apparent. In these calculations, we did not attempt to –nd the magnetic ground state, since the energy diÜerence between ferromagnetic and anti-ferromagnetic solutions is expected to be no more than a few meV.However, the HF calculations show that the antiferromagnetic solution is the most stable. Comparison of the spin densities yielded by the two theories (not shown here) shows them to be very similar, underlining the agreement as to the qualitative nature of the electronic ground state. This picture is in good accord with experiment: UPS measurements54,55 on the Ar-ion-bombarded (110) surface show that at low defect concentrations, a band of gap states exists, centred 2.3 eV above the valence band maximum. The gap state lies in the upper half of the gap, and moves up towards the bottom of the conduction band as the defect concentration increases.P.J. D. L indan et al. 143 Fig. 2 Densities-of-states, g(E) for spin-up (SC) and spin-down (SB) electrons, calculated using DFT, for the reduced 1]1 surface after structural relaxation.The energy zero is set at the topmost occupied state. In Table 2 we report the ionic relaxations. From DFT calculations, the stoichiometric results of Lindan et al.23 are reproduced for reference, along with the displacements on the reduced surface calculated without and with spin polarisation. To complete the table we show results for the reduced surface using unrestricted HF.10 Reduction of the surface produces very diÜerent relaxations compared with the stoichiometric case : the –ve-fold Ti now remains close to the bulk terminated position, the six-fold Ti relaxes into, rather than out of, the surface, and the in-plane oxygens move much further out of the surface.The results with spin polarisation are in good agreement for the two methods, and it is notable that the inclusion of spin polarisation makes little diÜerence to the relaxed geometry.An interesting feature is that the spin polarised DFT solution has small displacements of the in-plane oxygens along [001], thereby lowering the symmetry of the surface. This surprising result is a direct consequence of the nature of the excess-electron states, and is caused by the interaction of these states with the oxygen ions.The Ti(3d) orbitals rotate about [110] which, along with the oxygen displacements, reduces the overlap of the charge density on the oxygen ions with that of the Ti(3d) electrons. Table 2 Ionic displacements due to relaxation of the (110)1]1 surface (110)1]1 reduced (110)1]1 (a) (b) (c) (d) label [110] [16 10] [110] [16 10] [110] [16 10] [001] [110] [16 10] 1 0.09 0.00 [0.13 0.00 [0.11 0.00 0.00 [0.10 0.00 2 [0.12 0.00 0.00 0.00 0.01 0.00 0.00 [0.02 0.00 3 [0.09 0.00 » » » » » » » 4 0.11 [0.05 0.40 0.09 0.39 0.10 [0.04 0.39 » 5 0.11 0.05 0.40 [0.09 0.39 [0.10 0.04 0.39 » 6 [0.05 0.00 [0.02 0.00 [0.05 0.00 0.00 » » Labels refer to Fig. 1(c). The displacements are in and are from the bulk terminated positions.”, (a) stoichiometric surface ; (b) reduced surface without spin-polarisation ; (c) and (d) reduced spinpolarised surface. The results are from plane-wave pseudopotential DFT calculations (a)»(c) and HF calculations (d).144 Simulation of molecular processes on oxide surfaces Fig. 3 Plan view of the unrelaxed (100) surface. The 1]1 unit cell is represented by the rectangular box whilst the diamond box marks the c(2]2) supercell used.The c(2]2) adsorption geometry for the K ion is also shown. Non-irreducible atoms (those symmetry related to atoms in the asymmetric unit of the surface unit cell) are labelled with the suffix ìaœ. Perhaps the most striking feature of the DFT results is the energetics : for identical ionic con–gurations (the relaxed, stoichiometric positions minus the two bridging oxygen ions, a 16-ion system) we –nd that a spin-polarised solution with four more electrons spin-up than spin-down is 2.77 eV (1.1 J m~2) lower in energy than the spinpaired solution.After structural relaxation of both systems, the diÜerence in energy is 3.98 eV (1.6 J m~2).For the spin-polarised system the energy gained by relaxing the ions was 2.8 eV. IV. C Reduction by K adsorption on the (100) surface Controlled reduction of oxide surfaces may be achieved by adsorbing alkali-metal overlayers. Because of the intrinsic interest in the surface electronic structure and the technological importance of such overlayers a number of experimental studies have recently been performed on these systems.The adsorption of potassium on NiO56 and ZnO,57 caesium on NiO58 and sodium on MgO59 have been reported. The surface geometry of the adsorbed species is found to be strongly dependent on the nature of the substrate. Potassium on ZnO leads to an ordered p(2]2) overlayer whereas potassium or (00016 ) caesium on NiO and sodium on MgO forms islands around surface defects. On these metal oxides a strong interaction between the alkali metal and surface oxygen atoms is seen.Theoretical studies of these systems have been inhibited by the complex nature of the electronic structure. To our knowledge, no theoretical investigation of the nature of alkali-metal adsorption exists, although LMTO-DFT methods have recently been used60 to study transition metal overlayers on (110).TiO2 A number of recent experimental investigations have probed the geometry and electronic structure of potassium and sodium adsorption on surfaces.38,61h66 Low- TiO2 energy electron diÜraction (LEED) data indicate that K adsorption on the (100) surface at half-monolayer coverage results in an ordered c(2]2) geometry.61 This is con–rmed by surface-extended X-ray absorption –ne structure (SEXAFS) studies which also probe the local geometry of the K site.63 The electronic structure has been studied using photoelectron spectroscopy62 which revealed the existence of K-induced band gap statesP.J. D. L indan et al. 145 on titanium atoms some 1»2 eV above the valence band edge, as for the substoichiometric surface (see above).These states are apparently populated by charge transfered from the K overlayer to the substrate. In the present work, we have performed HF calculations on the adsorption of K on the (100) surface. The c(2]2) surface unit cell is shown in Fig. 3 and 4 in the TiO2 geometry obtained by relaxing the K position and top two layers of atoms , (O(1) O(2) , and The symmetry of the surface was retained ; symmetry equivalent atoms Ti(3) Ti(4)).are labelled with suffix a in Fig. 3. A number of adsorption sites for the K atom on the unrelaxed surface were studied. At the short bridge site (bridging and the O(1) O(2a)), ì 4-fold hollowœ site (equidistant from and and the atop site O(1), O(1a), O(2) O(2a)) (directly above one of the oxygen atoms), the energy was signi–cantly higher than that obtained for adsorption at the long bridge site (bridging Relaxation of the K O(1), O(2)).ion and the top two layers of the substrate results in the geometry displayed in Fig. 4. The optimal adsorption site is with the K atom 1.16 above the surface oxygen plane ” and equidistant (2.59 from the and atoms. This agrees well with the KwO ”) O(1a) O(2) bond length of 2.62 found in SEXAFS studies.63 The K is 0.06 nearer than ” ”Ti(4) The surface oxygen layer shows signi–cant relaxation as and move Ti(3).O(1) O(2a) closer together along the [001]-direction to increase the eÜective coordination of Ti(4) (see Fig. 4). The surface oxygen and the titanium layers (as de–ned by and Ti(3a) Ti(4)) move in opposite directions along the [010]-direction so that there is further increase in the eÜective oxygen coordination of the surface Ti sites.Similar relaxations have been observed on clean surfaces,17,21 the main diÜerence here being the tendency of TiO2 to have closer coordination to the K than Ti(4) Ti(3) . Some understanding of this rather curious adsorption site can be obtained by examining the charge and spin distributions at the surface which were estimated using a Mulliken population analysis,67 Table 3.The K atom is almost completely ionised, with 0.97 o e o being transferred to the surface. This charge is distributed over the (0.34 Ti(4) o e o) and the two surface O-sites (0.32 o e o). The site is almost unperturbed. The Ti(3) nature of this electronic state is even clearer in the spin populations.A single unpaired spin is generated on the site. The adsorption process can therefore be considered to Ti(4) involve charge transfer from the K atom to a speci–c surface Ti site generating a spin polarised Ti3` ion which, being negative compared to the lattice Ti4` ions, is less repulsive to the K. The density of states for this surface is compared with that in the bulk crystal and clean surface in Fig. 5. The broad valence band between [7.5 eV and the valence band maximum at 0.0 eV is dominated by orbitals and varies only in detail between the O2p bulk, clean surface and K-adsorbed surface. The clear signature of K adsorption is the band-gap state 1.4 eV above the valence band maximum, which is similar to the states Fig. 4 Relaxations of the substrate upon K adsorption. The K ion is shown in its optimised (”) position.146 Simulation of molecular processes on oxide surfaces Table 3 Mulliken population analysis of the charge (Qo e o ) and spin density ( o e o) of the c(2]2)- (100) K surface TiO2 atom Q Q[Qclean spin O(1) 9.49 0.32 0.00 O(2) 9.49 0.32 0.00 Ti(3) 19.34 0.06 0.01 Ti(4) 19.62 0.34 0.99 O(5) 9.48 0.02 [0.01 O(6) 9.48 0.02 0.00 O(7) 9.43 0.00 0.00 O(8) 9.44 0.02 [0.01 Ti(9) 19.19 0.00 0.01 Ti(10) 19.19 0.00 0.01 K 18.03 » 0.00 Fig. 5 HF densities of states projected onto Ti and O sites for (a) bulk and (b) the clean TiO2 , (100) and (c) (100)c(2]2) K surfaces. The shading indicates states derived from the induced Ti3` ion.P. J. D. L indan et al. 147 observed at 2»3 eV on the reduced (110) surface discussed above.As mentioned TiO2 above, this spin polarised state is almost entirely localised on the site. K-induced Ti(4) band-gap states are observed in photoelectron spectroscopy61h63 at 0.9 eV binding energy or 1»2 eV above the valence band maximum. V Water on (110) TiO2 Water interacts quite strongly with many oxide surfaces, and there have been many experiments on adsorbed water.In the case of MgO (001), there is now rather –rm evidence from experiment and theory that water is adsorbed in molecular form on the perfect surface, and that dissociation occurs only if defects are present.13 The situation on (110) is not so clear.68 Thermal desorption and UPS69 experiments have been TiO2 interpreted as showing that water is adsorbed dissociatively at sub-monolayer coverage.But very recently, high-resolution electron energy loss (HREELS) measurement70 have given evidence that dissociative adsorption occurs only at somewhat lower coverages, with molecular adsorption becoming dominant as monolayer coverage is approached. Recently we reported exploratory DFT-pseudopotential calculations of water adsorption, in which fully relaxed adsorption energies were obtained for a number of diÜerent con–guration.20 The simplest molecular candidate has the molecule coor- H2O dinated by its O atom to a –ve-fold surface Ti site ; the simplest dissociated con–guration has the water OH~ ion at the –ve-fold site (we call this a terminal hydroxy group), and the H] ion attached to a bridging O ion to form a second type of hydroxy ion (a bridging hydroxy group).With monolayer coverage and with the constraint of full symmetry (i.e. the molecular and OH~ axes point along the surface normal), a comparison of these two showed that the molecular con–guration is more stable, with Eads\ eV, compared with eV for the dissociated con–guration. However, the 0.82 Eads\0.45 symmetrical dissociated con–guration is only metastable, and by allowing the axes of the OwH groups to vary in a plane normal to the surface, distortion occurs, to produce hydrogen bonding between the two kinds of OH~ group; the resulting stabilisation increases to 1.08 eV.The conclusion is that for the con–gurations examined, disso- Eads ciation is energetically favoured. Support for dissociative adsorption at lower coverage comes from our recent dynamical DFT simulations22 in which the time evolution is obtained when a molecule coming from the gas phase reacts with the surface.We examined the adsorption of a single molecule on a surface doubled along [001] (i.e. with a 2]1 area), which corresponds to half-monolayer coverage. From several initial conditions studied we could not –nd a molecular adsorption site : once near one of the –ve-fold Ti sites, the molecule spontaneously dissociated on the surface, and the resulting con–guration was similar to that identi–ed in the monolayer-coverage static calculations.As an illustration of the time evolution, we show, in Fig. 6, three snapshots from one of the simulations, showing the spontaneous formation of the two types of surface hydroxy group mentioned above.After relaxation of the –nal dynamics con–guration, the distance between the O and H in the broken water bond was 1.8 . This is the typical length of a hydrogen bond and, ” indeed, there are clear signs of interaction between these atoms, as we describe shortly. One way to analyse the adsorbed species on a surface is to use an experimental technique which probes the vibrational spectra of those species.Henderson70 used HREELS techniques to do this for water adsorbed on (110). In interpreting these results, certain assumptions must be made as to the signals expected for the adsorbates. In particular, it is assumed that hydroxy groups contribute sharp, high-frequency peaks to the spectra.Thus, dissociation of water should be signalled by two such sharp peaks (the two hydroxy groups are in diÜerent environments) and an absence of the bond-bending frequency corresponding to the molecular vibrational mode in which the HwOwH angle varies.148 Simulation of molecular processes on oxide surfaces Fig. 6 Three snap-shots from a dynamical simulation of an molecule being dissociatively H2O adsorbed on the (110) surface.Large, medium and small spheres represent O, Ti and H TiO2 atoms respectively. Time proceeds from left to right, with an interval of ca. 0.2 ps between successive frames. Reproduced from ref. 22 with permission. We have performed equilibrium molecular dynamics (MD) simulations at low temperature, from which we have calculated the vibrational spectra of the hydrogen ions.Surprisingly, we –nd that only the terminal hydroxy group gives a well-de–ned high frequency signal. The bridging hydroxy group has a broad vibrational spectrum. These features are shown in Fig. 7. Direct examination of the vibrational motion shows that this broadening is due to the hydrogen-bonding interactions along the broken water bond: the thermal motion results in considerable variation in the length of this hydrogen bond, which consequently alters the vibrational frequency of the hydrogen ion.The high-frequency features we calculate appear to agree well with the HREELS spectrum in the corresponding frequency range. (We note that for technical reasons, these simulations employed an arti–cially large hydrogen mass of 3 u.This lowers the vibrational frequencies, but leaves the shape of the spectrum unchanged.) The conclusion from these calculations is that the vibrational signals of adsorbates may not be as anticipated from simple arguments. However, according to Henderson,70 a clear bond-bending signal is present at this coverage, a fact which our calculations do not explain. Of course, since we have only performed calculations at half-monolayer coverage in which a single molecule is treated explicitly, we cannot investigate intermolecular interactions which may be important at high coverage.Therefore, our most recent calculations aim to investigate the eÜect of intermolecular interactions on adsorption. We have taken the hydroxylated surface structure yielded by our previous MD studies of water adsorption22 and explored the adsorption of a second molecule on the surface.Again, we used several initial con–gurations, but the only adsorption site we found was over the second –ve-fold Ti site, i.e. the one not occupied by a hydroxyl group. At this site, adsorbs molecularly, and lies almost —at on the surface. This H2O con–guration, shown in Fig. 8 allows the formation of hydrogen bonds to both the surface bridging oxygen ion and the OH~ group which was produced by dissociation of the –rst molecule. Thus, from these calculations, we –nd that the most stable con–guration at monolayer coverage is one containing both dissociated and molecular water. Note that in this state, the separation of the molecule H and the terminal-hydroxy O, ca. 1.8 is very similar to the length of the hydrogen bond in the water dimer.72 In fact, ”, the geometry of the adsorbates is closely related to the dimer geometry. A metastable state, in which both molecules are intact, lies ca. 0.2 eV above this state. This remarkable result may provide an explanation of experiment, since it is clear that the adsorption state depends crucially on the interactions between molecules, and hence on the cover-P.J. D. L indan et al. 149 Fig. 7 Power spectra for the hydrogen ions in (a) the terminal hydroxy and (b) bridging hydroxy groups after water dissociation. The simulation was at 120 K. Note that the frequencies are shifted because of the use of a large hydrogen mass (see Section V). Note also the diÜerent y-scales used.Reproduced from ref. 22 with permission. age. In addition, the vibrational spectrum calculated from this surface contains a clear bond-bending signal, and is thus in good agreement with the HREELS data. The interaction of water with this surface is clearly a rather subtle matter, and the hydrogen bonds between water molecules and hydroxy groups are of crucial importance. Explicit examination of these interactions is vital if the correct adsorption state is to be determined.Also, since some of these states are of rather similar energy (within a few tenths of an eV), further re–nement and re-examination of our static calculations is warranted. We are, therefore, performing additional static-relaxation calculations using improved pseudopotentials. These pseudopotentials are of the ì ultrasoft œ type due to Vanderbilt.36 Their advantage is that smaller plane-wave cut-oÜs (i.e.smaller basis sets) can be used without sacri–cing accuracy, and this will allow us to go to larger supercells. We have used these pseudopotentials to investigate the wider range of adsorbate geometries shown in Fig. 9. In addition to the symmetric molecular (SM) and symmetric and unsymmetric dissociated (SD and UD) con–gurations studied before, we now include the unsymmetrical molecular con–gurations (UM1, UM2 and UM3) shown in the –gure.These new con–gurations were chosen on electrostatic grounds, which suggest that it should be favourable for the water O to coordinate to –ve-fold Ti and for both the H atoms to coordinate to O ions.This can be achieved in three ways: both Hs to the150 Simulation of molecular processes on oxide surfaces Fig. 8 Relaxed ionic con–guration for two water molecules adsorbed at monolayer coverage on the (110) surface. Ti and O are represented by light and dark large spheres, respectively, and the small spheres represent H. The relaxation was performed on the con–guration resulting from the dynamical simulation described in Section V.This is the lowest-energy con–guration found at monolayer coverage. (a) View looking down onto the surface. The 2]1 area of the simulation cell is outlined, and the dashed line indicates the boundary between the two 1]1 cells that form the surface. (b) View along Note the contraction along the H bond joining the molecule and the (116 0).OwH group. same bridging O (UM1); both Hs to two diÜerent in-plane Os (UM2); and one H to a bridging O and the other to the O of the water molecule in the neighbouring cell, which is one of its periodic images (UM3). We have calculated fully relaxed adsorption energies for these six con–gurations at monolayer coverage, using the three-layer slab.The values are systematically ca. 0.3 eV higher than those found in our previous static calculations, and we believe that part of the reason may be that the previous calculations were performed at a slightly larger lattice parameter. However, the adsorption energies for SM, SD and UD are in the same order as before, and the diÜerences between them are similar. The new molecular con- –gurations are all stable, with values of over 0.5 eV, and the most stable case, UM3, Eads has an which is only slightly less than that of the most stable dissociated con–gu- Eads ration.It is interesting to note that the molecular con–guration UM3 is very similar to the con–guration of the molecule in the mixed adsorption state found via the dynamics calculations at full coverage.P.J. D. L indan et al. 151 Fig. 9 Schematic representation of adsorbate geometries for water investigated using DFT staticrelaxation calculations. The labels are de–ned in the text. The 1]1 surface unit cell was used in all calculations, thus –xing the coverage at one monolayer. If one takes the static and dynamic calculations together, they suggest several things. First, dissociation is clearly favoured at half-monolayer coverage, but at full coverage there are signi–cant intermolecular interactions which make the molecular and dissociative adsorption energies similar.Another consequence of these interactions is that the molecule lies —at in the most stable molecular adsorption geometry, thereby forming a hydrogen bond along [001] with its neighbour.Our larger, two-molecule calculations allow additional stabilisation to occur through the contraction of this bond and, in fact, show that the lowest-energy geometry is the mixed dissociative/molecular state shown in Fig. 8. By analogy with the water dimer we may expect that this hydrogen bond contributes ca. 0.2 eV.71 Having demonstrated the importance of intermolecular interactions for water on this surface, we must note that the real adsorption geometry may have a periodicity other than the 2]1 case we have studied, but it is certain that it is not of 1]1 symmetry.Also, our results indicate that, at –nite temperatures, a given molecule will —uctuate between a dissociated and molecular state. Finally, we note that the lattice parameters of the substrate material may be crucial in determining the mode of adsorption.which also has the rutile structure, has larger a and c parameters and there- SnO2 , fore neighbouring molecules or OwH groups would be further apart than on the152 Simulation of molecular processes on oxide surfaces surface. The eÜect of this on the strength of the hydrogen bonds may result in very diÜerent behaviour. VI Conclusions Our ab initio calculations on have been used to study: the equilibrium structure of TiO2 the stoichiometric (110) surface ; the atomic and electronic structure of surfaces TiO2 reduced by removal of oxygen and by addition of potassium; and the molecular and dissociative adsorption of water at the (110) surface.We have shown that DFT and HF calculations give very similar predictions for the ionic displacements at the stoichiometric (110) surface, and that these are in semi-quantitative agreement with recent experimental measurements. The most important conclusion from our calculations on oxygen vacancies in both the bulk and at the surface is that spin unpairing is a dominant eÜect.Our DFT calculations show that the unpairing of spins leads to a substantial lowering of the total energy and to the appearance of the localised states in the band gap which are observed experimentally.Unrestricted HF calculations yield gap states with the same general characteristics. The two theoretical methods predict ionic displacements at the reduced (110) surface which agree very closely with each other. Our HF calculations show that K bonds to the (100) surface via charge transfer of an electron into a TiO2 particular surface Ti site.The resulting geometry is in excellent agreement with recent X-ray adsorption data. The gap states induced by surface reduction are similar to those observed on the non-stoichiometric surface. Finally, our static and dynamic DFT calculations on the adsorption of water on (110) are consistent with earlier –ndings that TiO2 dissociation gives stabilisation at half-monolayer coverage. However, our new dynamical calculations at monolayer coverage indicate a more complex picture in which water is present in both molecular and dissociated forms, which is consistent with recent HREELS measurements.The work of P.J.D.L. was done as part of a CCP3/CCP5 joint project supported by a grant from Cray Research Inc., and the work of S.P.B.and J.M. was supported by EPSRC grants GR/J34842 and 95563789 respectively. We also used local computer equipment provided by EPSRC grant GR/J36266. We are grateful for an allocation of time on the Cray T3D at EPCC provided by the High Performance Computing Initiative through the U.K. Car-Parrinello consortium and the Materials Chemistry consortium. Useful discussions with Prof.G. Thornton and Dr W. Temmerman are also acknowledged. This work has bene–ted from the collaboration within the Human Capital Mobility Network on ìAb initio (from electronic structure) calculation of complex processes in materialsœ. References 1 H. J. Freund, H. Kuhlenbeck and V. Staemmler, Rep.Prog. Phys., 1996, 59, 283. 2 G. Renaud, B. Villette, I. Vilfan and A. Bourret, Phys. Rev. L ett., 1994, 73, 1825. 3 G. Charlton, P. B. Howes, C. L. Nicklin, P. Steadman, J. S. G. Taylor, C. A. Muryn, S. P. Harte, J. Mercer, R. McGrath, D. Norman, T. S. Turner and G. Thornton, Phys. Rev. L ett., 1997, 78, 495. 4 M. J. Gillan, L. N. Kantorovich and P. J. D. Lindan. Curr. Opin.Solid State Mater. Sci., 1996, 1, 820. 5 M. Causa` , R. Dovesi, C. Pisani and C. Roetti, Surf. Sci., 1986, 175, 551. 6 C. A. Scamehorn, A. C. Hess and M. I. McCarthy, J. Chem. Phys., 1993, 99, 2786. 7 M. I. McCarthy, A. C. Hess, N. M. Harrison and V. R. Saunders, J. Chem. Phys., 1993, 98, 6387. 8 J. E. JaÜe, N. M. Harrison and A. C. Hess, Phys. Rev. B, 1994, 49, 11153. 9 C. Rehbein, N.M. Harrison and A. Wander, Phys. Rev. B, 1996, 54, 14066. 10 W. C. Mackrodt, E-A. Simson and N. M. Harrison, Surf. Sci., 1997, 384, 192. 11 C. A. Scamehorn, N. M. Harrison and M. I. McCarthy, J. Chem. Phys., 1994, 101, 1547. 12 L. N. Kantorovich, J. M. Holender and M. J. Gillan, Surf. Sci., 1995, 343, 221.P. J. D. L indan et al. 153 13 W. Langel and M. Parrinello, Phys.Rev. L ett., 1994, 73, 504. 14 I. Manassidis and M. J. Gillan, J. Am. Ceram. Soc., 1994, 77, 335. 15 M. J. Gillan, I. Manassidis and A. De Vita, Phil. Mag. B, 1994, 69, 879. 16 I. Manassidis, J. Goniakowski, L. N. Kantorovich and M. J. Gillan, Surf. Sci., 1995, 339, 258. 17 M. Ramamoorthy, D. Vanderbilt and R. D. King-Smith, Phys. Rev. B, 1994, 49, 16721. 18 M. Ramamoorthy, R. D.King-Smith and D. Vanderbilt, Phys. Rev. B, 1994, 49, 7709. 19 J. Goniakowski, J. M. Holender, L. N. Kantorovich, M. J. Gillan and J. A. White, Phys. Rev. B, 1996, 53, 957. 20 J. Goniakowski and M. J. Gillan, Surf. Sci., 1996, 350, 145. 21 P. J. D. Lindan, N. M. Harrison, J. Holender, M. J. Gillan and M. C. Payne, Surf. Sci., 1996, 364, 431. 22 P. J. D. Lindan, N. M. Harrison, J.M. Holender and M. J. Gillan, Chem. Phys. L ett., 1996, 261, 246. 23 P. J. D. Lindan, N. M. Harrison, J. A. White and M. J. Gillan, Phys. Rev. B, 1997, 55, 15919. 24 M. C. Payne, M. P. Teter, D. C. Allan, T. A. Arias and J. D. Joannopoulos, Rev. Mod. Phys., 1992, 64, 1045. 25 R. Dovesi, C. Pisani and C. Roetti, Hartree-Fock Ab Initio T reatment of Crystalline Systems, Volume 48, Springer-Verlag, Berlin, 1988. 26 R. Jones and O. Gunnarsson, Rev. Mod. Phys., 1989, 61, 689. 27 J. P. Perdew in Electronic Structure of Solids œ91, ed. P. Ziesche and H. Eschrig, Akademie Verlag, Berlin, 1991. 28 J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, D. J. Singh and C. Fiolhais, Phys. Rev. B, 1992, 46, 6671. 29 V. Heine, Solid State Phys., 1970, 24, 1. 30 Cohen and V.Heine, Solid State Phys., 1970, 24, 37. 31 A. Zupan and M. Causa` , J. Quantum Chem., 1995, 56, 337. 32 M. D. Towler, A. Zupan and M. Causa` , Comput. Phys. Commun., 1996, 98, 181. 33 P. Reinhardt, B. A. Hess and M. Causa` , J. Quantum Chem., 1996, 58, 297. 34 L. J. Clarke, I. and M. C. Payne, Comput. Phys. Comm., 1992, 72, 14. Sã tich 35 G. Kresse and J. Furthmué ller, Phys. Rev.B, 1996, 54, 11169. 36 D. Vanderbilt, Phys. Rev. B, 1990, 41, 7892. 37 R. Dovesi, V. R. Saunders, C. Roetti, M. Causa` , N. M. Harrison, R. Orlando and E. Apra` , CRY ST AL 95 Userœs Manual, University of Torino, 1996. 38 V. E. Henrich and P. A. Cox, T he Surface Science of Metal Oxides, Cambridge University Press, Cambridge, 1994. 39 V. E. Henrich and R. L. Kurtz, Phys. Rev. B, 1981, 23, 6280. 40 Z. Zhang, S.P. Jeng and V. E. Henrich, Phys. Rev. B, 1991, 43, 12004. 41 C. H. Wang and Y. S. Xu, Surf. Sci., 1989, 219, L 537. 42 M. Tsukada, H. Adachi and C. Satoko, Surf. Sci., 1983, 14, 113. 43 S. Munnix and M. Schmeits, J. V ac. Sci. T echnol. A, 1987, 5, 910. 44 E. Pickett, Rev. Mod. Phys., 1989, 61, 433. 45 J. Zaanen and G. Sawatzky, Prog. T heor. Phys. Suppl., 1990, 101, 231. 46 A. Svane and O. Gunnarsson, Phys. Rev. L ett., 1990, 65, 1148. 47 Z. Szotek, W. M, Temmerman and H. Winter, Phys. Rev. B, 1993, 47, 4029. 48 M. D. Towler, N. L. Allan, N. M. Harrison, V. R. Saunders, W. C. Mackrodt and E. Apra` , Phys. Rev. B, 1994, 50, 5041. 49 N. M. Harrison, in Computer Modelling in Inorganic Chemistry, ed. C. R. A. Catlow, Academic Press, San Diego, 1997. 50 N. M. Harrison, B. G. Searle and E. A. Seddon, Chem. Phys. L ett., 1997, 266, 501. 51 J. B. Goodenough, Prog. Solid State Chem., 1972, 5, 145. 52 A. M. Ferrari and G. Pacchioni, J. Phys. Chem., 1995, 99, 17010. 53 P. A. Cox, T ransition Metal Oxides, Oxford University Press, Oxford, 1992. 54 V. E. Henrich, G. Dresselhaus and H. J. Zeiger, Phys. Rev. L ett., 1976, 36, 1335. 55 H. R. Sadeghi and V. E. Henrich, J. Catal., 1988, 109, 1. 56 S. Kennou, M. Kamaratos and C. A Papageorgopoulos, V acuum, 1990, 41, 22. 57 D. Purdie, C. A. Muryn, S. Crook, P. L. Wincott, G. Thornton and D. A. Fischer, Surf. Sci. L ett., 1993, 290, L680. 58 S. Kennou, M. Kamaratos and C. A Papageorgopoulos, Surf. Sci., 1991, 256, 312. 59 H. Onishi, C. Egawa, T. Aruga and Y. Iwasawa, Surf. Sci., 1987, 191, 479. 60 L. Thien-Nga and A. T. Paxton, Phys. Rev. B, submitted. 61 R. Casanova, K. Prabhakaran and G. Thornton, J. Phys. Condens. Matter, 1991, 3, S91. 62 P. J. Hardman, R. Casanova, K. Prabhakaran, C. A. Muryn, P. L. Wincott, G. Thornton, Surf. Sci., 1991, 269/270, 677. 63 K. Prabhakaran, D. Purdie, R. Casanova, C. A. Muryn, P. J. Hardman, P. L. Wincott and G. Thornton, Phys. Rev. B, 1992, 45, 6969. 64 P. W. Murray, N. G. Condon and G. Thornton, Surf. Sci., 1994, 323, L281.154 Simulation of molecular processes on oxide surfaces 65 A. G. Thomas, P. J. Hardman, C. A. Muryn, H. S. Dhariwal, A. F. Prime, G. Thornton, E. Roman and J. L. Desegovia, J. Chem. Soc., Faraday T rans., 1995, 91, 3569. 66 R. Heise and R. Courths, Surf. Rev. L ett., 1995, 2, 147. 67 R. S. Mulliken, J. Chem. Phys., 1955, 23, 1833. 68 V. E. Henrich and P. A. Cox, T he Surface Science of Metal Oxides, Cambridge University Press, Cambridge, 1994. 69 M. B. Hugenschmidt, L. Gamble and C. T. Campbell, Surf. Sci., 1994, 302, 329. 70 M. A. Henderson, Surf. Sci., 1996, 355, 151. 71 S. S. Xantheas and T. H. Dunning, J. Chem. Phys., 1993, 99, 8774. Paper 7/02103H; Received 26th March, 1997