Faraday Discuss., 1997, 106, 1»40 Introductory Lecture Computer modelling as a technique in solid state chemistry C. Richard A. Catlow, Lutz Ackermann, Robert G. Bell, Furio Cora` , David H. Gay, Martin A. Nygren, J. Carlos Pereira, German Sastre, Ben Slater and Phillip E. Sinclair T he Royal Institution of Great Britain, 21 Albemarle Street, L ondon, UK W 1X 4BS We describe the role of computer modelling techniques in solving problems relating to both the structural and dynamical properties and to the reactivity of solids.Recent illustrative examples are taken from the –elds of microporous materials, ionic and transition metal oxides and molecular crystals. 1 Introduction Atomistic computer modelling techniques now play a major role in both physical and biological sciences. In this paper, we aim to show how simulation and modelling methodologies may be used as tools to solve complex problems in solid state chemistry in a way that is complementary to experimental investigations.We illustrate some of the current possibilities in this –eld by recent application from our laboratory, which will describe the modelling of structures and bonding in inorganic materials, of surfaces and morphologies of both oxide and molecular materials, of defect structures of solids and of solid state reactivity with a strong emphasis on catalytic applications.We will also show how modelling is beginning to have an impact on our understanding of the hydrothermal synthesis of silicates»traditionally a very empirical –eld. We aim to exemplify by these applications, together with the many others presented in this volume, the expanding range and scope of this –eld. 1.1 Methodologies Methodological aspects of the –eld are described in several recent reviews1h5 and in many of the articles within this volume; indeed, the methods employed are increasingly the standard techniques of computational chemistry and physics adapted to the particular problems posed by solids.We therefore present only a brief summary. Methods based on interatomic potentials have a wide range of applicability in solid state chemistry. Energy minimisation (EM) methods (often referred to, in the solid state context, as static lattice techniques) have been used for many years, particularly in modelling crystal structures6 and defect energies and structures ;7 in the latter case, methods based on the procedure originally pioneered by Mott and Littleton8 have enjoyed particular success.Despite the simplicity and limitations of minimisation methods, they are robust and widely applicable. The techniques may be extended to include free energy minimisation by including the calculation of vibrational entropies using lattice dynamical theory.7,9 Such methods have been fruitfully applied to studying phase transitions, especially in silicate minerals.9 12 Computer modelling as a technique in solid state chemistry Molecular dynamics (MD) methods have been extensively employed for modelling diÜusion in solids where atomic/ionic transport is rapid, as in fast ion conductors,10 although for slower diÜusion, which is more typical of solids, such methods are inappropriate as the number and range of particle displacements is insufficient in the timescale of an MD simulation, which in current work rarely exceeds 1 ns ; for slower diÜusion we must use methods based on transition state theory which identify saddle points for atomic jump mechanisms. Application to molecular diÜusion within microporous solids has been a particularly fruitful –eld for MD, as will be discussed later in this article, although again the techniques are limited to rapidly diÜusing species, with approaches based on transition state theory again being used for slower diÜusions. MD has also been used to probe rotational dynamics in solids,11 and in a limited number of cases to model directly phase transitions.12 Monte Carlo (MC) methods have proved to be of particular value in studying sorption in solids.13 Grand canonical Monte Carlo (GCMC) methods are being increasingly used in this –eld, although simpler methods which blend MC, MD and EM techniques have proved eÜective in identifying low energy sites for sorption in microporous solids.14,15 A diÜerent and important range of applications relate to kinetic MC studies, in which the MC ìmovesœ are atomic jumps whose frequencies have been calculated using static lattice methods.Several applications have been reported by Murch16,17 to alloy diÜusion; this approach was also used in modelling oxygen diÜusion in the widely used solid state electrolyte, yttrium doped More recently, a successful applica- CeO2 .18 tion of this approach to the diÜusion of benzene in zeolite Y was reported by Auerbach et al.19 The nature of the interatomic potentials employed in such calculations has been extensively reviewed elsewhere,20,21 and further details and discussion will be given in later sections of this article and in other papers in the volume.The type of potential model employed depends on the nature of the bonding in the solid.Molecular mechanics potentials, whose conceptual base is the covalently bonded solid, and which contain terms depending on the deviation of bond lengths and angles from equilibrium values, is obviously appropriate for covalent systems. Representations of polarisability should, wherever possible, be included. The shell model22 has provided a simple and robust way of modelling such eÜects in oxide and halide systems, although the limitations of these potentials, particularly in calculating phase transition energies, has become apparent. 23,24 Parametrisation of potentials can be undertaken using empirical approaches which –t variable parameters to observed properties of model compounds. Alternatively (and increasingly) parameters are calculated directly from electronic structure studies of clusters or periodic arrays of atoms, examples of which are given in ref. 25»27. As the –eld of computational solid state chemistry moves towards increasingly predictive applications, there is a growing need for increased accuracy in interatomic potential models. Explicit electronic structure techniques have of course been used since the early days of solid state theory.The major development in recent years has been the growth, due to the development in both methodology and hardware, of the applicability of high quality ab initio methods, although there is a continuing role for semi-empirical techniques. Local density functional (LDF) methods are being widely used for periodic systems (especially the plane-wave pseudopotential approaches of the type pioneered by Car and Parrinello28), and contemporary work increasingly emphasises the use of gradient corrected density functionals rather than those based on the simpler local density approximation.Hartree»Fock (HF) methods are also extensively used, where the periodic boundary conditions methods, available in the CRYSTAL code,29 have had notable impact.At the time of writing, calculations using ab initio methods on unit cells of ca. 50 atoms are becoming increasingly practicable (although larger systems such as the zeolite ZSM-530 with 192 atoms in the unit cell have been studied). Cluster methods are often more appropriate for modelling defects, although these may be described using periodicC. R.A. Catlow et al. 3 supercells. When clusters are employed, care must be taken, whenever possible, to embed them in a suitable representation of the surrounding lattice. Interfacing of the cluster with its embedding matrix may pose serious problems in such calculations. Recent work in the –eld is described in their proceedings by Sierka and Sauer31 and Sherwood et al.32 With the continuing growth of the applicability of electronic structure techniques, can we see them as replacing interatomic potential based methods? In some cases, this will unquestionably be the case.An example is the case of defect structures in relatively simple ceramic oxides, e.g., MgO, and which could only be modelled by TiO2 Al2O3 , Mott and Littleton8 type methods in the early 1980s, but which are now within the range of good quality electronic structure calculations. However, the more general answer is that there will be a continuing role for interatomic potential based methods as the –eld moves to more complex systems. Moreover, there are many problems, e.g., the modelling of the structures of complex silicates, and of diÜusion in microporous materials, where such methods are the most accurate and appropriate.A –nal comment is needed on hardware, the growth of which constantly expands the horizons of the –eld, and where todayœs workstations are more powerful than yesterdayœs supercomputers. It is hard to predict how the –eld will develop. The growth of usable massively parallel systems is certainly one major trend, which computational solid state chemistry is indeed increasingly exploiting.There is no doubt, however, that the exponential growth in computer power will continue for the foreseeable future. In the sections which follow, we will highlight recent applications in which both interatomic potential and electronic structure methods have been used. 2 Modelling structures Modelling of the structure at the atomic level of a solid is of course a necessary –rst step in any solid state simulation study. If we are unable to reproduce the structure with acceptable accuracy, then we can have little con–dence in subsequent applications to the study of, e.g., lattice dynamical or defect properties.Moreover, when modelling techniques are applied to surfaces or to the structure of glassy phases, it is essential to ensure that they adequately reproduce the structures of the corresponding bulk crystalline phases.Additional incentives for accurate structure modelling are provided by the need for guidance from computationally based methods to the solution of the structures of complex crystals, e.g., microporous materials. Indeed, computational methods are now becoming powerful structure solving tools in the crystallography of complex inorganic materials,5 as shown by the work of Gorman et al.,33 described in this volume.As has been described elsewhere,34 application of modelling methods to crystal structures has passed through three distinct phases: in the –rst, the methods are used to reproduce known structures ; in the second, we re–ne approximate structures ; while in the third (and most ambitious class of applications) we aim to predict structures from little initial knowledge of the arrangements of the atoms within the unit cell. The techniques employed in the –rst two categories are energy minimisation, starting from a known approximate structure ; several versatile computer codes, WMIN,35 METAPOCS, 36 THBREL37 and GULP38 are available for such calculations. Good examples of structure reproduction are found in the work of Bell et al.39 on the monoclinic phase of the microporous zeolite silicalite and the recent detailed exploration of the structure of high silica zeolites by Henson et al.40 who achieved a good measure of agreement between calculated and experimental studies.The work of Shannon et al.41 provides an excellent illustration of the role of modelling in re–ning approximate structures. In this impressive study, modelling techniques were used to generate a structure for the high silica zeolite Nu 87 from an approximate initial structure (which was unable to be re–ned adequately from the available high resolution powder diÜraction data). The –nal,4 Computer modelling as a technique in solid state chemistry energy minimised structure, however, led to a fully satisfactory re–nement of the crystallographic data.Structure prediction is essentially an exercise in global optimisation. The most widely used procedures involve simulated annealing which uses a Monte Carlo algorithm to explore the potential-energy surface of a system and to identify low energy regions, which are then re–ned by EM techniques. Often the simulated annealing calculation uses a simpli–ed cost function to assess the suitability of the structure without recourse to a full lattice energy calculation.An alternative approach uses genetic algorithms which permute (breed) structures in such a way that they evolve towards low cost functions. Recent work, which solved the crystal structure of is a good illustra- Li3RuO4,42 tion of the latter approach.Examples of the use of simulated annealing techniques will be given in the paper of Gorman et al.33 in this volume. A challenging new problem in the –eld of structure modelling is provided by the mesoporous silica structures which will now be described. 2.1 Computer modelling of mesoporous materials Since the discovery of the –rst ordered mesoporous molecular sieves with pore diameters in the range 25»100 there has been a rapid growth in research into these fasci- ”,43,44 nating materials. Much eÜort is currently focused on two areas : (1) rational design of structured inorganic mesopores, with the aim of understanding synthetic mechanisms and of exploiting this knowledge to obtain materials with a desired composition, structure and properties, and (2) the modi–cation of these materials to introduce catalytically active centres into the mesopore structure.Examples of the latter include the substitution of transition metal atoms, such as Ti, into the pore framework of siliceous MCM- 41.45,46 In many cases the short-range structure of these mesopores is not well understood.In particular, the original MCM-41 phases, prepared by liquid crystal templating, form highly regular hexagonal arrays of pores, and were shown to have pore walls as thin as two T atoms in width. However, there was no evidence for any crystallographic order, other than that caused by the two-dimensional ordering of the mesoporous channels, and it has been generally accepted that the pore walls consist largely of amorphous silica, a conclusion which is supported by the presence of a high concentration of terminal silanol groups.It is nevertheless reasonable to suppose that there exists some degree of ordering of tetrahedra. In our initial approach to the problem of meso- SiO4 pore structure, we have therefore followed two approaches: –rst we have designed regular zeolitic mesopore model structures in which, as in conventional zeolites, each T atom is connected to four others by an oxygen bridge.Secondly, we have used a completely amorphous silica surface, still containing four-coordinate T atoms, but replete with terminal SiOH groups. Each model has been subjected to energy minimisation calculations, both as the purely siliceous isomorph, and with titanium substituted into the pore framework. The calculations employed are the standard lattice and defect energy minimisation available in the GULP38 program using the potential parameters reported by Jentys and Catlow.47 2.1.1 Results and Discussion.As –rst postulated by Smith and Dytrych,48 it is possible to construct an in–nite series of hexagonal zeolite structures with uniform channels of 6n T atoms, where n\2 represents the AlPO-5 structure and, as it transpired, n\3 is equivalent to the topology of VPI-5.We have extended this series to obtain structures with eÜective pore diameters of greater than 25 and have so far examined 30-, 36- and ”, 42-ring structures (see Fig. 1). Each framework was constructed from tetrahedra SiO4 and subjected to full energy minimisation using the program GULP.38 The principal features of the optimised structures are summarised in Table 1.C. R. A. Catlow et al. 5 Fig. 1 Energy minimised hypothetical 36-ring hexagonal structure The lattice energies of the hypothetical frameworks, given in Table 2, would seem not to be in—uenced by the density of the material.This result is at variance with trends observed in microporous zeolites where there is a correlation between framework density and lattice energy per unit, and suggests a levelling oÜ in lattice energy TO2 when all, or nearly all, of the tetrahedra lie on the internal surface of the SiO4 materials»a –nding that is consistent with experimental work of Navrotsky et al.49 The incorporation of titanium into the framework was also investigated by substituting a single titanium atom onto a T site in the 36-ring structure and carrying out an energy minimisation of the surrounding lattice.The energy of substitution fell within the range calculated by Jentys and Catlow47 for the various T sites of silicate, in fact within 1 kJ mol~1 of that of the most stable site.The four TiwO bond lengths were in the range 1.74»1.80 ”. Table 1 Summary of energy minimised ìzeolitic mesoporeœ structures T atoms in ring a/” c/” free diameter/” space group 30 27.96 8.36 21.0 P6 6 2m 36 32.40 8.51 25.3 P6/mcc 42 37.03 8.36 30.0 P6 6 2m Table 2 Lattice enthalpies of ì zeolitic mesoporesœ relative to other silica polymorphs material lattice energy/kJ mol~1 per TO2 quartz 0 silicalite 11 30-ring 24 36-ring 27 42-ring 266 Computer modelling as a technique in solid state chemistry Fig. 2 Energy minimised amorphous silica surface, showing the substituted group TiO3OH Amorphous pore walls have also been modelled by taking a dense silica glass structure, obtained by simulated annealing, cleaving out a thin layer, approximately three T atoms in thickness, and satisfying all dangling SiwO bonds with hydroxy groups.Part of an energy minimised structure is shown in Fig. 2 where a Ti atom has again been substituted into the pore wall. In this case the TiwO bond lengths varied from 1.77 to 1.92 (although, because of the presence of the OH group, the substitution energies ” cannot be directly compared with those in the zeolite frameworks). 3 Characterisation of bonding in solids Theoretical investigations aimed at understanding the basic interactions in the solid state chemistry of complex materials, like transition metal oxides, have until recently relied on simpli–ed and parametrised schemes, often based on extended-Hué ckel or tightbinding Hamiltonians.50,51 With the improvement in hardware and software tools now available, ab initio quantum mechanical studies are becoming increasingly feasible, and can be employed to examine in detail the properties of bonding in the solid state, in the absence of empirical parameters and simpli–ed theoretical frameworks.As an example of the kind of insight and detailed information now possible, we report in this section on Hartree»Fock (HF) and density functional (DF) calculations on three binary oxides with stoichiometry i.e., and In both MO3 , MoO3, WO3 ReO3 .schemes, the model adopted to represent the bulk materials is based on periodic boundary conditions. In the HF calculations, the electronic distribution of the solid is described in terms of crystalline orbitals, obtained as a linear combination of localised functions, or atomic orbitals (AOs), associated with lattice sites.The HF code employed is CRYSTAL.52 DF calculations are based on the full-potential linear muffin-tin orbital (FP-LMTO) code developed by Methfessel et al.,53 within the local density approximation (LDA) and the Hedin»Lundqvist exchange-correlation potential ; basis functions are in this case Hankel functions augmented by numerical solutions of the radial Schro é dinger equation in the muffin-tin spheres.In the following discussion we focus on several distinct features of the solid state chemistry of transition metal oxides. 3.1 Nature of the transition metal The three oxides examined have common structural features, all based on the ReO3 structure, which is formed by a network of corner-sharing octahedra, with all the MO6C.R. A. Catlow et al. 7 oxygen ions in a two-coordinate, bridging position between adjacent metal sites. To illustrate the nature of the bonding in these systems, we make reference –rst to the highest symmetry, cubic phase. In Fig. 3 we report the band structure of project- ReO3 , ed onto its components in the Re and O basis sets, as obtained from the DF calculation. 54 The valence band (VB) is formed mainly by O 2p states, although they are mixed with the Re 5d orbitals. The conduction band (CB) has a dominant contribution from the Re 5d orbitals, while the Re 5d orbitals lie at a higher energy; this energy t2g eg dependence is typical of d AOs in an octahedral –eld.It is important to note the presence of three levels along the C»X direction, denoted in Fig. 3 with the labels ìAœ and ìBœ (in the CB) and ì1œ (in the VB), which are pure Re pure Re and pure O 2p t2g, eg levels respectively, with no mixing of orbitals. These levels are —at in k-space, and have been previously denoted as superdegenerate.55 Covalent interactions between the relevant O 2p and Re 5d AOs are symmetry forbidden in the cubic phase. Levels A, B and 1 are therefore non-bonding between metal and oxygens, which interact in the solid via their ionic charge.However, the importance of covalence in the bonding is shown by the bottom levels of the VB, which have a high k-dispersion and contribution from both Re and O. The hybridisation pattern of metal and oxygen levels can be interpreted, in classical terms, as a back-donation of electrons from the –lled O2~ AOs to the empty 5d AOs on the metal.The solution for and is very similar to that described MoO3 WO3 here for ReO3 . The energy gap between the superdegenerate levels A and 1, *E(A[1), is related to the value of the crystal –eld in the metal and oxygen positions : the decrease of the eÜective nuclear charge (taking into account the number of core electrons) of the transition metal ion from Re to Mo and W aÜects the relative energy of levels A and 1.The energy of level A, which as shown in Fig. 3 is a metal d AO, is more sensitive to the Fig. 3 Band structure of cubic projected onto its atomic contributions, along the C»X»M ReO3 , direction of reciprocal space.The bands labelled ìAœ and ìBœ (in the conduction band) and ì1œ (in the valence band) are pure Re pure Re and pure O 2p states respectively. t2g, eg8 Computer modelling as a technique in solid state chemistry nuclear charge of M than the energy of level 1, which is an O 2p state. *E(A[1) is therefore higher in and than it is in Moreover, Mo and W have a MoO3 WO3 ReO3 .similar eÜective nuclear charge, but the orbitals involved are the 4d in Mo and the 5d in W; the gap *E(A[1) is lower in than it is in This feature is important, MoO3 WO3 . because smaller values of the gap *E(A[1) correspond to a higher degree of covalence in the bonding.55,56 Hence, *E(A[1) to an extent measures the tendency towards covalent interactions : the lower *E(A[1), the more covalent the solid.The trend in the balance between covalent and ionic bonding is con–rmed by the population analysis data, as obtained from the HF study, reported in Table 3. A higher value of the bond population corresponds to a higher degree of covalence in the bonding, as is also (qb) con–rmed by the lower values of the net ionic charges on metal and oxygens which accompany the higher values. qb 3.2 Number of valence electrons We now examine the symmetry reductions from the cubic phase: experimental evidence shows that the cubic structure of is stable at all temperatures, while and ReO3 MoO3 are stable in lower symmetry phases, in which the transition metal ion lies oÜ- WO3 centre in its coordination octahedron.For simplicity we focus here on the cubic]tetra- Table 3 Equilibrium bond distances and Mulliken population analysis within each octa- MO6 hedron in bulk and as obtained from the HF study.(The indexing of atoms ReO3 , MoO3 WO3 , is the same as adopted in the text). bond distance bond population atom multiplicity net charge r(MwO)/” qb(MwO) cubic ReO3 Re ]3.31 O 6 [1.10 1.85 ]0.03 cubic b-MoO3 Mo ]2.99 O 6 [1.00 1.89 ]0.01 cubic WO3 W ]4.09 O 6 [1.36 1.88 [0.03 tetragonal b-MoO3 Mo ]2.74 O(1) 1 [0.70 1.70 ]0.17 O(2) 4 [1.02 1.87 [0.05 O@(1) 1 [0.70 2.30 ]0.01 SOTa 6 [0.92 2.86 ]0.03 tetragonal WO3 W ]3.79 O(1) 1 [0.97 1.63 ]0.21 O(2) 4 [1.41 1.90 [0.01 O@(1) 1 [0.97 2.37 [0.00 SOTa 6 [1.26 1.93 ]0.03 orthorhombic a-MoO3 Mo ]2.43 O(1) 1 [0.44 1.64 ]0.29 O(2) 1 [0.75 1.69 ]0.10 O(3) 2 [1.25 1.94 [0.02 O@(2) 1 [0.75 2.21 [0.02 O@(3) 1 [1.25 2.24 [0.00 SOTa 6 [0.81 1.94 ]0.06 a Average M»O value in an octahedron. MO6C.R. A. Catlow et al. 9 gonal distortion, which corresponds to a displacement of the transition metal ion along the [001] crystallographic direction, towards one of its six nearest neighbour oxygen ligands. Fig. 4 shows the band structure of the tetragonal phase of after a ferroelectric ReO3 , movement of the Re sublattice by 3% of the lattice parameter along the [001] direction. The mixing of the superdegenerate levels A and 1 in the distorted structure is no longer symmetry forbidden, and it gives rise to a p bonding level in the VB and to a p antibonding level in the CB.A similar but smaller change also aÜects the levels of r symmetry.In the distortion, therefore, a new covalent bond is formed between transition metal and oxygen, which involves AOs of both r and p symmetry along the MwO directions. The relative stability of the cubic and tetragonal phases depends on the number of valence electrons (n) : in and only the bonding level 1 in the VB is –lled, MoO3 WO3 , and the distortion is energetically stable ; in there is, in addition, one antibonding ReO3 electron in the CB, which opposes the distortion.It has been shown54 that the electronic degrees of freedom in the two phases are isoenergetic for a number of electrons in the CB, When the cubic phase is stable and vice versa for n0\0.98 o e o. n[n0 0\n\n0 where the distortion is favoured, in agreement with experiment.With reference to the population analysis of Table 3, we see that the increase in covalence is restricted to the interaction between the transition metal and only one of its six nearest neighbour oxygens, the one towards which the metal displaces [labelled O(1) in Table 3], while the interaction with the four equatorial oxygens [O(2) in Table 3] is mostly unchanged with respect to the cubic phase.The above feature allowed us to label Fig. 4 Projected band structure of after the Re displacement along the [001] direction. ReO3 Bands ìAœ and ì1œ are the same as in Fig. 3; now they have a completely diÜerent atomic contribution, showing the eÜect of the newly formed p M 2p bond in the distorted structure. t2gwO10 Computer modelling as a technique in solid state chemistry the short MwO pairs in and as tungstyl and molybdenyl groups.57 Fur- MoO3 WO3 thermore, the displacement of the transition metal leaves a longer MwO bond, and creates a more layered bulk structure.57 Mechanical properties can also be shown to depend on n:54 in Table 4 we report the calculated bulk moduli (B) for the oxides examined.The very high B of re—ects the ReO3 population of the RewO p antibonding orbital in the CB: by increasing the applied pressure, the stress acting along the RewO directions is very high.In practice, the solid reacts to this stress by diverting the applied pressure from the RewO towards the OwO direction ; this eÜect will cause a rotation of the oxygen octahedra, and lead to a better exploitation of the empty spaces of the corner-sharing network, in agreement with the experimentally observed pressure-induced phase transition.In and the MoO3 WO3 , eÜect of pressure is absorbed by a change in the interlayer distance ;57 in both cases the eÜect of the structural change is such that the value of B decreases. 3.3 Structural requirements So far, we have examined the in—uence of the transition metal ion on the overall properties of the materials ; this is not the only relevant aspect in describing their solid state chemistry : several transition metal oxides can in fact exist in diÜerent polymorphs, in which the metal ion is the same, but the medium-range arrangement of the ions is diÜerent.To investigate this feature we examine the two polymorphs of the MoO3: (the same structure as discussed in Sections 3.1 and 3.2), and the ReO3-like b-MoO3 layered and thermodynamically stable a-MoO3 .In the crystal structure of (Fig. 5) we identify a direction (a), perpendicular a-MoO3 to the layers, along which the octahedral units are disconnected from one MoO6 another; within the layers, the connectivity of the octahedra occurs in one direc- MoO6 tion (b) by common edges, so as to form zigzag rows, and in the perpendicular (c) by common corners only.Because of this anisotropy, oxygens in have a diÜerent a-MoO3 local environment: the outermost oxygen in each layer has only one nearest molybdenum atom; where the octahedra are connected by corner-sharing, the oxygens are two-coordinated, as in while oxygens along shared edges are three- a-MoO3 , coordinated.In the following discussion, we refer to n-coordinated oxygens as O(n). With reference to the population analysis data reported in Table 3, we clearly see that oxygens with a diÜerent coordination number in the solid behave like diÜerent chemical species.58 In particular, the Coulomb –eld created by only one nearest Mo in the O(1) positions is not strong enough to stabilise an O ion with the formal charge of [2.As a consequence, the interaction between Mo and O(1) is much more covalent than that between MowO(2) and MowO(3); the particularly low value of the net charge on O(1) Table 4 Calculated bulk modulus (B/GPa) for the optimised structures of and ReO3 , MoO3 WO3 , as obtained from our DF and HF study phase DF HF HF]corra cubic ReO3 304 282 320 cubic b-MoO3 » 279 316 cubic WO3 254 257 281 tetragonal WO3 » 137 174 a Results obtained estimating the contribution of electron correlation a posteriori, based on a functional of the HF equilibrium density.C.R. A. Catlow et al. 11 Fig. 5 Network of edge and corner sharing octahedra in showing the layering of MoO6 a-MoO3 , the structure along the a crystallographic direction con–rms this feature.The short MowO(1) pair in behaves again as a single a-MoO3 unit, or molybdenyl group. Another important structural feature to examine in is the connectivity a-MoO3 between layers. In the inter-layer region, we –nd molybdenyl groups facing one another, and connecting in a zip-like manner (see Fig. 6). The low net charge on O(1) contributes to decreasing the repulsion between layers, and stabilises the layered polymorph with respect to As shown above, is more ionic than the net charge on b-MoO3.WO3 MoO3 ; O(1) in a layered polymorph of analogous to would therefore be higher, WO3 a-MoO3 hence destabilising this structure.58 3.4 Long-range Coulombic forces Finally, we show the eÜect of long-range forces, by comparing ferro- and antiferroelectric (FE and AFE respectively) displacements of W in Both structures have WO3 .the same short-range chemical features and medium-range structural requirements, as commented in Sections 3.1»3.3 above, but diÜer in the long-range order in the solid. The movement of one ion in a solid is equivalent to creating a net dipole in the structure ; in the FE phase all dipoles are aligned, while in the AFE the dipoles in neighbouring cells are anti-aligned. Only the Coulombic, long-range interaction will diÜerentiate the two structures.In Fig. 7 we report the energy change for FE and AFE displacements of W in we see there that short- and long-range interactions play a comparable role in WO3; determining the solid state chemistry of WO3 .From all the previous arguments we must conclude that the bonding in the solid state has several components. We may distinguish between a chemical (short-range) contribution, which depends on the nature of the transition metal and on the number of valence electrons, a structural (medium-range) contribution depending on the connectivity of the octahedra in the solid, and a Coulombic (long-range) contribution MO6 which takes into account the ionic interaction.The equilibrium con–guration is due to a12 Computer modelling as a technique in solid state chemistry Fig. 6 DiÜerence electron density map (solid minus isolated Mo6` and O2~ ions), in the plane through the outermost molybdenyl groups [MowO(1)] in adjacent layers of high- a-MoO3 , lighting the zip-like connection of layers and the short O(1)wO(1) equilibrium distance.The zigzag dashed line in the centre of the plot serves as a guideline to separate the molybdenyl groups belonging to each of the two facing layers. The lines correspond to positive (»»), negative (» » ») and zero (» … » …) diÜerence. The interval between the isodensity lines is 0.005 e a0~3.subtle balance of all these forces. Accurate quantum mechanical calculations, of the kind described here, can help us in –nding their relative importance in the structures under investigation, and we can use them to rationalise experimental data on these complex materials. Fig. 7 Dependence of the total energy of tetragonal on ferro- and antiferro-electric WO3 (L) (Ö) displacements of W; the oxygen sublattice is kept –xed in the cubic phase.represents the ZW displacement of W (in fractional coordinates) from the centre of the octahedron. WO6C. R. A. Catlow et al. 13 4 Modelling of the surfaces of complex materials Surfaces play a central role in the chemistry and physics of solids. The details of the surface structure control many of the observed properties of the system, such as growth/ dissolution, catalysis and resistivity as in gas sensors.To understand these phenomena we require knowledge of the surface structure, and the surface structure in diÜerent environments. Modelling plays a complementary role to experiment in this –eld, due to the complexity of making accurate measurements of the surface structure, and difficulty in preparing clean surfaces. 4.1 Techniques The technology of modelling surfaces is very similar to that used in modelling bulk systems, but with diÜerent boundary conditions. The same potential models used for the bulk can generally be used for the surface although results for surfaces tend to be more dependent on the quality of the potential than for those of bulk systems, owing to the reduction of symmetry at the surface.The new boundary conditions are vacuum or solvent over a surface with 2D (two-dimensional) periodic boundaries parallel to the surface normal. These boundary conditions require the use of more parameters to describe the surface uniquely. Specifying the surface normal does not identify a unique crystal face. A position along the surface normal, or shift, is needed as well.The normal and shift are sufficient for surfaces that are planar or smooth cuts. Molecular crystals (to be described later) require special care since molecules should not be fragmented to create a surface. In addition to the structure we can compute the surface energy and the attachment energy. The former is the energy required to cleave the crystal to create the surface.The latter is the energy released when a growth slice is attached to the surface. These quantities can be used to predict the crystal morphology.59 Use of the surface energy implies thermodynamic control of the morphology. The conditions under which use of attachment energy is appropriate are ill de–ned, but they are generally more suited to morphologies controlled by kinetic factors. 4.2 Strategies for modelling surfaces There are three ways to model a 2D system. One is a 3D solution that uses a 3D array of slabs and gaps. The second is to compute the total energy using 2D periodic boundary conditions. The third uses a 3D lattice and computes the energy of the atoms while keeping track of the side of a cutting plane on which they are located.This procedure does not permit the relaxation of the surface and therefore is of limited applicability. While this partitioning technique can be used for quick calculations of the attachment energy, it is not of general applicability for modelling surfaces and hence will not be discussed further. Many simulations, both atomistic and quantum mechanical, are carried out using a 3D array of slabs. The total energy of the system is the energy of the atoms in the cell interacting with their lattice images.However, there are restrictions on when a 3D slabgap approach can be used. The simulation cell must not have a dipole moment normal to the slab or have one created during the simulation, as the dipole moment would interact across the gap distorting the adjacent surface.To guarantee that this problem does not arise, a slab must have a plane of mirror symmetry located in the middle of the slab. The other restriction is that the top and bottom of the slab should be equivalent, otherwise it is not certain to which face the computed properties apply. The other full 2D approach as used in the program MARVIN59 (and in the earlier MIDAS code60) is more general and can be used even when the restrictions on the 3Dz y x region 1 region 2 14 Computer modelling as a technique in solid state chemistry slab-gap apply. The energy of the cell is computed with a 2D lattice summation (with an adaptation of the Ewald sum to 2D).In the 2D approach, the bulk of the crystal can be treated as an array of –xed atoms.This two-region approach is illustrated in Fig. 8. For energy minimisation or dynamics, the region 2 atoms are –xed to reproduce the force exerted by the bulk crystal. The region 1 atoms are allowed to move either to minimise the total energy or dynamics can be introduced using the MD technique. The full 2D, two-region technique also allows for docking calculations. These permit the attachment or incorporation of impurities or molecules at or near the surface.The eÜects of impurities on crystal growth or other surface properties can be determined from this type of calculation. 4.3 Types of surface We can divide surfaces into two classes : polar and non-polar. Non-polar faces are the easiest and most reliable to model, since no corrections need to be made to compute the surface energy or attachment energy.Accurate modelling of polar faces requires special attention, as it is necessary to neutralise the dipole in some way, either mechanical, environmental or chemical. Mechanical neutralisation involves moving some ions from one face to another to counteract the dipole ; chemical neutralisation removes some of the charge by attaching new charged species, such as hydroxy groups, to the surface ; while environmental neutralisation is eÜected by polarisation of the surrounding solvent.We can also distinguish between three diÜerent types of face, depending upon the structure of the unit cell in the direction of the surface normal (Fig. 9). Type 1 faces are non-polar and type 3 faces are polar, regardless of the position of the cutting plane.Type 2 faces are either polar or non-polar depending upon where the cutting plane is. Crystals composed of molecules present special problems to modelling of surfaces. A planar cut through the bulk crystal to create the surface usually passes through molecules located near the surface. Using this plane will leave molecular fragments on the Fig. 8 The two-region surface simulation cell.Region 1 is relaxed to minimise the total energy while region 2 is –xed and represents the bulk crystal.C. R. A. Catlow et al. 15 Fig. 9 The three types of surface. Type 1 has no dipole moment normal to the surface, wherever it is cleaved. Type 2 surfaceœs dipole moment depends upon the cleavage plane, some of which have no dipoles. Type 3 faces have a dipole, regardless of where the cleavage plane is located.surface, which is not chemically or physically realistic. Creation of such fragments is avoided in the MARVIN code. We now give three examples of recent applications of surface modelling. The –rst relates to the key problems of surface hydration of an inorganic material. The second describes recent quantum mechanical studies of transition metal oxide surfaces.The third describes recent surface simulations of molecular materials. 4.4 Hydroxylation of the corundum basal plane surface : modelling of surface chemical reactions Surfaces commonly undergo stabilising chemical reactions. The surface ions and molecules are exposed to the environment and a clean surface will contain coordinatively unsaturated ions and dangling bonds.For surfaces having a net dipole, the dipole must be cancelled to create an electrostatically stable structure, as discussed above. A chemical reaction that reduces the surface charge may stabilise a cut that is unstable as a clean surface. In this section we discuss an important and topical example, namely the hydroxylation of the basal plane surface of corundum. Replacing an oxygen with a hydroxide ion halves the charge; so hydroxylating an oxygen terminated cut reduces the charge by a factor two which, depending on the interplanar spacing, may cancel the net dipole.61 For the oxygen terminated surface this is the case and other examples Al2O3(0001) where hydroxylation cancels a net dipole are oxygen terminated (111) surfaces of rocksalt structured oxides or (001) surfaces of perovskite oxides.The clean, aluminium terminated basal plane surface has been studied by Al2O3 several authors using both atomistic potentials59,62 and periodic Hartree»Fock methods.63 Both types of calculation predict that the top layers collapse inwards in a dramatic fashion. These large relaxations have recently been observed in experiments by Guenard et al.using X-ray scattering at grazing incidence.64 The aluminium terminated surface exposes three-coordinate aluminium ions and half of the aluminium sites are empty, as shown in Fig. 10. A surface exposing both dangling bonds and ions with unsaturated coordination can be expected to be highly reactive. Replacing the topmost oxygen layer with hydroxy ions for the oxygen terminated surface, not only cancels the dipole but gives rise to a surface without any dangling bonds and no low coordinated ions, as shown in Fig. 11. In sharp contrast to the aluminium terminated surface, the hydroxylated surface shows hardly any relaxation, as illustrated in Fig. 12. We now explore the energetics of creating hydroxylated surfaces. For a clean surface, the surface energy only includes the work needed to cleave the perfect crystal.When the16 Computer modelling as a technique in solid state chemistry Fig. 10 Unrelaxed and relaxed aluminium terminated surface of Half of the aluminium Al2O3 . sites are empty and upon relaxation the topmost aluminium relaxes inwards. Small light spheres represent aluminium ions ; larger, darker spheres represent oxygen ions surface undergoes a chemical reaction, the reaction enthalpy must also be included into the surface energy.The ideal way of calculating such energies would of course be to use quantum mechanical methods rather than interatomic potential based procedures. The computational cost of quantum mechanical methods often rules out this possibility and we are faced with the problem of incorporating reaction enthalpies into calculations using interatomic potentials.In the case of hydroxylation, we need to calculate the energy for creation of the surface together with the reaction : Osurf 2~ ]H2O]2OHsurf ~ Such energies of reaction are normally calculated by introducing a hypothetical gasphase reaction, which should then be chosen so that the heat of reaction can be esti- Fig. 11 The hydroxylated surface of where all oxygens have been replaced by hydroxides Al2O3 , leaving a surface with no low coordinated ions or dangling bonds. There is almost no surface relaxation.C. R. A. Catlow et al. 17 Fig. 12 Changes in interplanar spacings upon relaxation. The middle column shows the inter- (”) planar spacing for the unrelaxed surface.The left hand column is the relaxed aluminium terminated surface and the right hand column the hydroxylated surface. mated for all steps, either from calculations or thermochemical data. A Born»Haber cycle commonly used for investigating hydroxylation is therefore Osurf 2~ ]Ogas 2~ H2Ogas]Ogas 2~]2OHgas ~ 2OHgas ~ ]2OHsurf ~ In this reaction we have the problem that O2~ does not exist in the gas phase, which means that we have to de–ne the second electron affinity of oxygen.There have been several previous estimates of this quantity using thermochemical data; for example, a value of 8.27 eV was estimated65 which when used in calculating the surface energy for the hydroxylated surface gives a negative value of [1.12 J m~2.66 A negative surface energy for a clean surface would mean that it would be energetically favourable to maximise the surface area.It is more plausible that negative energies could occur when the surface undergoes a chemical reaction, for which it is also possible that further reaction might be kinetically inhibited. Nevertheless such values must be viewed with considerable scepticism. An alternative procedure is to calculate the enthalpy for the hydroxylation reaction using the experimental heat of formation for the corresponding metal hydroxide67 where this is available.We wish, however, to develop a reliable and general method of calculating the surface energy for reactive surfaces using a procedure that does not rely on hypothetical species, such as gaseous O2~. We will now outline such a scheme together with recent results for the case of the hydroxylated surface.We –rst note that for Al2O3 reactions containing only neutral species it is much easier to calculate the reaction enthalpy using quantum mechanical methods than it is for reactions involving negatively18 Computer modelling as a technique in solid state chemistry charged ions.We therefore choose a gas-phase reaction that is charge neutral. It must also be possible to calculate the enthalpies for surface to gas and gas to surface processes. Starting from the oxygen terminated surface we move half of the oxygens to the other side of the crystal neutralising the dipole. Reacting the resulting surface with water gives rise to the hydroxylated surface.To calculate the surface energy of the hydroxylated surface we begin by calculating the surface energy for this half-layer oxygen terminated surface. Starting from this surface and, keeping the above considerations about the gas-phase reaction in mind, a suitable Born»Haber cycle for calculating the reaction enthalpy is given in Table 5 together with the calculated energies for each reaction step.In the –rst step, we pull oÜ a growth layer from the oxygen terminated surface. To have well de–ned energetics, this layer must have no dipole moment across it ; hence the choice of as the repeat unit. In the next step we separate this layer into clusters. Al4O6 After reacting the clusters with water we put them back together again into a hydroxylated layer.Finally we drop this layer back onto the surface. The enthalpy for reacting the cluster with water has been calculated using gradient corrected density func- Al4O6 tional techniques.68 All other steps have been calculated using a consistently derived set of atomistic potentials.66 If we normalise the total energy given in Table 5 per unit surface area, we obtain a value of [9.95 J m~2.On adding this value to the surface energy for the half-terminated oxygen surface (13.40 J m~2) we calculate a value of 3.45 J m~2. Relaxing the hydroxylated surface we obtain for the –nal surface energy a value of 3.40 J m~2. Comparing this value with that of [1.12 J m~2 obtained using the cycle employing gaseous O2~ illustrates the care that is needed in constructing the appropriate thermochemical cycle.We also note that the relaxed surface energy for the aluminium terminated surface is 2.52 J m~2. We should stress that there still remain uncertainties in our results arising from the chosen interatomic potential parameters and from the geometries used in the DFT calculation on the gas-phase reaction. Nevertheless, the present calculation suggests that creation of the unhydroxylated surface costs less energy than of the hydroxylated.However, our earlier calculations suggest that, once created, the hydroxylated surface will be stable ; indeed we calculated the water desorption energy to be at least 5 eV. There is no inconsistency or paradox here. Removal of water from the hydroxylated surface exposes the high energy oxygen terminated rather than the lower energy aluminium terminated surface and hence is highly endothermic.It would be of great interest to prepare the hydroxylated surface and to examine its structure. 4.5 Quantum mechanical calculations on surfaces : bonding in surface layers Employing similar arguments to those used in Section 3 to examine bulk materials, we can try to understand the eÜect of surfaces on the bonding properties of transition metal Table 5 Calculated energy changes (*E) for the component steps in the Born»Haber cycle for the hydroxylation of the surface of Al2O3 *E/eV Al4O6 (surface)]Al4O6 (layer) 60.97 Al4O6 (layer)]Al4O6 (cluster) 19.05 Al4O6 (cluster)]3H2O]Al2O3Al2(OH)6 (cluster) [13.85 Al2O3Al2(OH)6 (cluster)]Al2O3Al2(OH)6 (layer) [22.62 Al2O3Al2(OH)6 (layer)]Al2O3Al2(OH)6 (surface) [69.33 total energy: [25.78C. R.A. Catlow et al. 19 oxides, a factor of great importance in improving our understanding of the catalytic properties of these materials. In this section we focus on perfect and sur- WO3 MoO3 faces, which we have represented by a slab model, a 2D array of atoms of –nite thickness, described using the same computational conditions employed for the HF study of the respective bulk systems.As we have shown in Section 3, both tetragonal and the two polymorphs of WO3 can be described as layer structures. The layering orientation in bulk provides a MoO3 preferential direction of cleavage for the surfaces : the S001T for the materials, ReO3-like and the S100T for Let us examine the structure –rst : stable surfaces a-MoO3 .ReO3 must be charge neutral and must not have net dipoles perpendicular to the surface.69 In the structure each oxygen is two-coordinated and shared between two transition ReO3 metal ions ; if we imagine the surface as created by cleaving the bulk material along a predetermined plane, the requirements summarised above mean that one half of the oxygens belonging to the cleavage plane must be left in each of the two facing surfaces created.In practice this causes a doubling of the unit cell, and a )2])2 reconstruction of the surface with respect to the bulk structure, as observed experimentally.70 The S001T surfaces of and of are illustrated in Fig. 13(a). The transition metal WO3 b-MoO3 ions on the surface are alternately –ve- and six-coordinated ; in the former case only the longest MwO bond in the bulk arrangement has been cleaved on the surface, thus representing a relatively small perturbation to the bulk arrangement.No major structural and electronic rearrangement occurs on the surface, and the overall bonding properties remain the same as in the bulk materials ; the only relevant diÜerence is the Fig. 13 Schematic representation of the S001T surface of and (a), and of the S102T b-MoO3 WO3 surface of (b). Atoms are labelled here as in the text. WO320 Computer modelling as a technique in solid state chemistry shortening of the surface molybdenyl and tungstyl groups and the corresponding increase in their covalence. Results of the population analysis for the surfaces examined are reported in Table 6.[Note that in Table 6 and in the following discussion we have used the following notation for labelling surface atoms: whenever an index was used in the bulk structure (see Table 3), this has been retained in the surface notation, to highlight the original position of surface species in the bulk; an index ì s œ has been added to the ions directly exposed onto the surface.Metal atoms are further labelled with index 5 or 6, according to their coordination number on the surface.] The –ve-fold coordinated Mo(5s) and W(5s) lie below the plane of the four equatorial oxygens [O(2s)], by which they are eÜectively screened from a direct exposure onto the surface. Similarly, the S100T surface of is cleaved along the direction of layering of a-MoO3 the bulk material, and its bonding remains unchanged.These surfaces are the stable ones for the oxides examined; in cleaving the solid along any other plane we will cut stronger MwO bonds, which will cause a stronger perturbation, and major electronic and structural rearrangements on the surface. To exemplify this feature we examine in more detail the S10nT surfaces of tetragonal these correspond to S001T platelets, WO3: separated by kinks every n unit cells.In particular, we examine one member of this family, the S102T surface [Fig. 13(b)]. For the same requirements of surface stability expressed above, W ions on kinks are four-coordinated (following the notation introduced above, we label atoms on kinks with index ìkœ).Three of their four O ligands are in turn two-coordinated, and continue into the S001T platelets ; they are accordingly labelled O(ls) and O(2s). The fourth oxygen is instead coordinated only to W(k), and we refer to it as O(2k). As we might expect, most of the rearrangement involves W(k) and O(2k): in the unrelaxed structure, O(2k) corresponds to one of the ionically bound equa- Table 6 Equilibrium bond distances and Mulliken population analysis for the surfaces examined of and MoO3 WO3 bond distance bond population atoma multiplicity net charge r(MwO)/” qb(MwO) a-MoO3 S100T Mo(s) ]2.40 O(1s) 1 [0.39 1.65 ]0.31 O(2s) 1 [0.76 1.69 ]0.09 O(3s) 2 [1.25 1.94 [0.02 b-MoO3 S001T Mo(6s) ]2.64 O(1s) 1 [0.44 1.70 ]0.36 O(2s) 4 [1.04 1.86 [0.09 Mo(5s) ]2.67 O@(1s) 1 [0.69 1.70 ]0.18 O(2s) 4 [1.04 1.87 [0.05 WO3 S001T W(6s) ]3.65 O(1s) 1 [0.66 1.59 ]0.35 O(2s) 4 [1.44 1.91 [0.03 W(5s) ]3.73 O@(1s) 1 [0.96 1.63 ]0.20 O(2s) 4 [1.44 1.90 [0.02 WO3 S102T W(k) ]3.31 O(2k) 1 [0.79 1.65 ]0.27 O(1s) 1 [1.10 1.60 ]0.06 O(2s) 2 [1.46 1.87 ]0.01 a Labelling of surface atoms is explained in the text.C.R. A. Catlow et al. 21 torial oxygens, O(2) (compare with Table 3).The relaxation on the surface involves the displacement of O(2k) along the diagonal of the three W(k)wO(s) directions ; the metal ion on the kink, W(k), is therefore in a distorted tetrahedral environment. The energy gained in the relaxation of O(2k) is 2.55 eV, con–rming that the electronic rearrangement accompanying the geometric relaxation is substantial : as we could expect from the analysis of the bulk system, the now singly coordinated O(2k) left free on the surface behaves in a completely diÜerent way from both the corresponding bulk species and from two-coordinated O(2).In particular, covalence in the bond between W(k) and O(2k) is very pronounced (as con–rmed by the bond population reported in Table 6) and grows partially at the expense of the W(k)wO(1s) interaction, which is still covalent, although to a much lesser extent than in the bulk material.The properties of bonding close to the kink are therefore completely diÜerent from those of the bulk species. This –nding is of general validity : the stronger the perturbation created by the surface, the higher the electronic rearrangement and the structural reconstruction on the surface. We expect the chemical reactivity of the exposed surface species to be greatly in—uenced by the above features.Finally, in Fig. 14 we show the electrostatic potential created by the surfaces described above. This is a most important property of surfaces, which controls the adsorption of gas-phase molecules. We stress an important common feature : a negative –eld prevails only above O(1s) species, while the –eld in the region above O(2s) and O(3s) ions is positive. Although the latter have the highest negative charge among the oxygens, the –eld which they create is compensated and dominated by the high positive charge of the neighbouring cations.In adsorption processes, the positive end of molecules would therefore be attracted by the O(1s) end of surface tungstyl and molybdenyl groups, which is very important in determining the chemical properties of the exposed surfaces, for example in catalysis and in gas-sensing devices.The –eld in the region next to kinks is much more intense and structured, and would therefore drive polar or polarisable molecules preferentially towards these surface defects. 4.6 Simulation of the surface of molecular crystals Previous sections have summarised the rapid progress in the –eld of computational surface science, particularly in the area of ionic surfaces and related surface mediated processes. However, in the case of molecular crystals, computational expense precludes all but the simplest molecular crystals from treatment via periodic all-electron or semiempirical methods.71,72 Therefore, the majority of molecular surface or periodic slab simulations reported in the literature are based on interatomic potential techniques.Molecular crystals are distinct from other classes of crystals in that their stability is controlled primarily by dispersive forces, and their packing arrangement by short-range repulsion. This precludes these crystals from treatment via periodic Hartree»Fock methods (except in the case of strongly hydrogen bonded solids : see, for example, ref. 71 and 72) and plane wave based techniques. The majority of electronic molecular surface studies reported are based upon cluster techniques,72 where short-range and polarisation eÜects dominate. We have previously shown how ionic surfaces may reconstruct in order to eliminate permanent perpendicular surface dipoles, giving rise to a divergent surface potential.In molecular crystals, unlike ionic or semi-ionic materials, the molecular fragments often exhibit permanent intrinsic dipoles, giving rise to a surface dipole, but the methods by which the dipole may be eliminated for these materials are less well understood.Molecular surfaces may contain highly polarisable functional groups, which are clearly aÜected by solvents and may undergo substantial relaxation. In order to simulate the polar morphology of urea, Docherty et al.72 reported a study where the potentials and charge22 Computer modelling as a technique in solid state chemistry Fig. 14 Maps of the electrostatic potential generated by three of the surfaces examined: a-MoO3 S100T (a) ; S001T (b) and S102T adjacent to the kink (c).The electrostatic potential b-MoO3 WO3 map for S001T is not reported, due to its close resemblance to plot (b), relating to WO3 b-MoO3 S001T. Positive (»»), negative (» » ») and zero (» … » …) potentials are shown; the interval between isopotential lines is 0.005 au (1 au\27.21 eV).distribution were scaled, based on PHF and cluster studies, in order to distinguish between polar mirror planes. However, their technique does not include surface relaxation, yet the issues of relaxation and potential model modi–cation are interdependent. Relaxation of molecular surfaces is characterised by a more shallow potential-energy surface minimum in contrast to electrostatically dominated ionic surface relaxation.Two classical based codes are commonly employed to explore molecular surface processes : MARVIN (discussed previously59) and HABIT.73 These programs diÜer in their treatment of surfaces, in that HABIT does not consider the eÜects of surface relaxation, which has been shown to be of importance by George et al.74 in the case of urea. In molecular crystals the dispersive and electrostatic forces are required to be computed particularly accurately : MARVIN employs a full 2D Ewald75 sum, an approach that is critical when considering polar molecules.C.R. A. Catlow et al. 23 Several previous studies of the surfaces of molecular crystals have employed the periodic bond chain technique: Berovitch-Yellin et al.76,77 have successfully shown how to manipulate additives to in—uence crystal morphology and Docherty et al.72 have reported a study of the eÜects of solvents on morphology.In order to explore these factors and kinetic eÜects quantitatively and comprehensively, we have developed two extensions to the functionality of MARVIN. MARVIN has previously been used successfully to predict the morphology and surface properties of a number of organic and inorganic crystals59,74,78,79 and has been extended to incorporate automated docking procedures via MC and MD techniques. The addition of these methods represents a signi–cant development in the capabilities of surface simulation.The automated docking scheme allows one to explore efficiently the conformational space of single or multiple additive species. This is of great value in the study of the eÜect of additive/impurities on the surface/morphology.The velocity Verlet algorithm80,81 and Gear82,83 predictor-corrector schemes have been implemented, allowing one to study, for example, kinetic eÜects, surface diÜusion, crystal growth and dynamic interfacial activity. Previously, surface MD studies have generally been carried out using a periodic slab approach (see, for example, Boek et al.84), where, as noted earlier, the simulation box has to be constructed such that a mirror plane is not present, as this may give rise to a dipole perpendicular to the surface. One can combine the MD and MC techniques to enable efficient exploration of conformational space, by using low energy con–gurations generated from MC simulations as initial conformers in MD simulations.We have exploited both these techniques to investigate factors aÜecting the crystal growth of the important molecular crystal aspartame. Aspartame is widely used within the food industry as an arti–cial sweetener, but extraction of robust crystals can be difficult and expensive. Recently, we have been addressing this problem, and carried out simulations using the techniques outlined earlier in this review.Initially, an energy minimisation of the bulk aspartame crystal was carried out, using the CFF91 force –eld,85 extracted from the DISCOVER module supplied by MSI.86 To eliminate unnecessary computational eÜort, we selected only the most signi–cant torsional terms from the force –eld.(Full details of the potentials used and further details regarding the simulations reported here will be given in subsequent publications.) The lattice parameters were reproduced within 1% of those reported in the crystal structure, 87 an excellent agreement for a complex molecular crystal. This relaxed cell was used to generate the four most morphologically signi–cant surfaces : subsequent relaxation of these surfaces showed minimal changes in surface energies and attachment energies, indicating a predicted morphology qualitatively similar to that produced by a Donnay»Harker88 analysis.As in the case of urea, aspartame has polar faces which appear in the morphology. To address the eÜect of polar solvents, on apolar and polar surfaces, we have simulated saturated aspartame, with an overlayer of water.One feature of molecular surfaces that complicates this kind of treatment is the problem of de–ning the solvent surface-accessible volume. This problem is compounded, in the case of aspartame, where the polar surfaces are highly porous and the apolar surfaces are corrugated. Having de–ned the accessible surface, relaxed bulk water was grafted onto the surface.The bulk water was obtained by a 500 ps MD run on a cell containing 1024 atoms at 298 K using an SPC89 potential model. The resulting cell box was relaxed at constant pressure using the static lattice method available in the GULP38 code. Solvent molecules with anomalously small overlaps with the solute of \1.2 ” were eliminated. The resulting surface and interface was relaxed using the MARVIN code, and this structure used as the initial con–guration in a surface MD run.Fig. 15 shows a snapshot taken after 2 ps, including 1 ps equilibration. After 5 ps, the water is seen to diÜuse in well de–ned channels into the aspartame surface. Longer simulations24 Computer modelling as a technique in solid state chemistry Fig. 15 Snapshot of the aspartame(110)/water interface after 2 ps of MD also show dynamic exchange of crystal water with surface water. We are currently carrying out extended MD runs to establish quantitatively the in—uence of the solvent on surface and attachment energies. In summary, the advent of full 2D Ewald surface MD allows one to probe dynamical surface phenomena of crystal surfaces that are intractable or approximate by other methods.This capability has general implications and in this article we have highlighted its particular relevance to simulating solvent eÜects in a complex molecular crystal. 5 Simulation of defects There is a long history in the application of defect modelling techniques in solid state chemistry and physics.5,7 Indeed, the extensive range of studies published in the 1970s and 1980s established Mott»Littleton related techniques as quantitative simulation tools in the –eld.Defect calculations continue to be a fertile and productive application area, as we will illustrate by recent applications90,91 concerning the mechanisms of dissolution of water in upper mantle minerals»a key problem in current mineralogy, asC.R. A. Catlow et al. 25 Table 7 Calculated energies of water dissolution in Mg2SiO4 energy/eV phase reaction (I) reaction (II) olivine 3.71 0.46 b phase 1.18 [1.10 the upper mantle is known to contain substantial quantities of dissolved water.92,93 The broader questions raised by the mechanisms of solution of water in oxides and silicates have been discussed in ref. 94 and 95. The upper mantle minerals are magnesium silicates of general formula and have the olivine structure at lower (Mg0.88Fe0.12)SiO4 temperatures, and a spinel like structure at higher pressures and an intermediate b phase at intermediate pressures. The transition zone in the upper mantle is known to be associated with a phase transformation from the olivine to the b phase.Wright and Catlow90,91 have clearly shown that the dissolution of water in both the olivine and b phase may be eÜected by two mechanisms: MgMg]2OO]H2O][VMg(OH2)]]ìMgOœ (I) 2FeMg~]2OO]H2O]2FeMg]2OH~]12 O2 (II) in which we use the Kroé ger»Vink notation and where the –rst reaction involves the formation of a defect cluster containing a magnesium vacancy with two protonated neighbouring atoms.The second reaction involves reduction of iron(III) (substituting for magnesium) to iron(II), with charge compensation produced by the protonation of neighbouring oxygen ions. The calculated energies for these reactions in both olivine and the b phase are given in Table 7. We note that the energetics of the iron(III)/iron(II) redox reaction is low in both phases, and the calculated energies of water dissolution are lower in the b phase compared with olivine»a result that is in accord with experiment. There is little doubt that the mechanisms shown in reactions (I) and (II) are the main processes responsible for dissolution of water in the upper mantle minerals. 6 Modelling of diÜusion As noted in the Introduction, MD techniques have been widely used in simulating diÜusion in solids and a particularly fruitful –eld of application has concerned the transport of sorbed molecules, especially hydrocarbons within microporous hosts.Many chemical reactions undergone by hydrocarbons when using zeolites as catalysts are very diÜerent when compared to the same reactions carried out in the gas phase.96h98 The hydrocarbons have, of course, to diÜuse through these structures in order to reach the active site.The large number of zeotypes99 allows us, in principle, to select among a wide variety of structures in order to favour the yield of a particular reaction. Indeed, shape selectivity100 is commonly a diÜusion controlled process which in—uences product distribution in reactions such as aromatic isomerisation among others.The understanding of more recently investigated phenomena like secondary shape selectivity and inverse shape selectivity among others illustrates the importance of hydrocarbon diÜusion in in—uencing reactivity in zeolites.101,102 DiÜusion measurements of hydrocarbons in zeolites have been made for several decades,103 although often many experimental variables other than the microscopic diffusion itself in—uence the measured diÜusion coefficients. In addition to direct microscopic experimental techniques like PFG-NMR104 and quasi-elastic neutron scattering (QENS),105 computer simulations can give very good estimates of self-diÜusivities.10626 Computer modelling as a technique in solid state chemistry In this section, we address the study of ortho- and para-xylene diÜusion in CIT-1,107 the –rst microporous material with 10 MR (membered rings) and 12 MR channels.DiÜusion of these two C8 aromatics has been widely studied in 10 MR zeolites like ZSM-5 where the size of the channels precludes the diÜusivity of the ortho isomer.108 DiÜusion of these isomers in a structure also containing 12 MR, through which both isomers can diÜuse, brings new features, in particular the competition of both isomers for the diÜusion path through the 12 MR channels when only one of them can diÜuse through the narrower 10 MR channels.We have, therefore, undertaken atomistic MD simulations to model the diÜusion of ortho- and para-xylene in purely siliceous CIT-1. The calculations require the speci–- cation of a potential-energy function which allows the calculation of its –rst derivatives from which the forces acting on each atom can be obtained.The force –elds employed include four terms as follows : Vtotal\Vzeolite]Vxylene]Vxylenehxylene]Vzeolitehxylene where the zeolite potential includes Buckingham, three-body and Coulombic terms; the intramolecular xylene potential includes two-body, three-body, four-body and Coulombic terms; and the xylene»xylene and zeolite»xylene potentials have Lennard-Jones and Coulombic terms.A more detailed description of the potentials can be found in ref. 109. All the MD simulations were performed using the general purpose DLñPOLY2.0110 code in its parallel version running on a 512 PE CRAY-T3D using 128 processors. The simulations proceed by assigning initial velocities to the atoms and then solving Newtonœs equations of motion using a –nite time step by means of the standard Verlet algorithm.80 Timesteps of 1 fs and an equilibration temperature of 500 K were employed in the simulations.The system comprises a 4]2]4 macrocell of purely silica CIT-1 and a maximum of 32 xylene molecules with a total of 3264 atoms to which periodic boundary conditions are applied throughout the MD simulation.Simulations of 100 ps were carried out within the NVE ensemble at the three diÜerent xylene loadings of 32 ortho-xylene, 8 ortho-xylene and 8 para-xylene molecules in the CIT-1 macrocell. Additionally, activation energies for the diÜusion path of each isomer in each channel were calculated by means of the solids diÜusion path facility implemented in the MSI Catalysis 6.0 software package.111 This simple but eÜective procedure drags a molecule through a channel in a microporous structure, calculating the interaction energy of the molecule with its host.The results of the three simulations yield three basis pieces of information, namely activation energies, the diÜusion path followed by the molecules through the zeolite channels, and self-diÜusion coefficients.The activation energies are shown in Fig. 16. It can be seen that ortho- and para-xylene can diÜuse without important restrictions through the 12 MR channels as shown by the low activation energies (6.21 and 7.03 kcal mol~1 for the para and ortho isomers respectively). Large interactions may be calculated when ortho-xylene diÜuses through the 12 MR channels due to an unoptimised orientation of the molecule in the channel (see Fig. 16) but a reorientation of the methyl groups towards the channel direction produces a low energy diÜusion path.This feature has important implications as to the possibility of rotation of the isomers, a point to which we return later. In the 10 MR channels a diÜerent result arises, as the activation energy for the para isomer of 20.93 kcal mol~1 is still consistent with diÜusion, whereas the 110.31 kcal mol~1 for the ortho isomer make its diÜusion impossible through the narrower channel.It is interesting to note that some penetration of the ortho isomer in the 10 MR channels is allowed at the intersection between the 10 MR and the 12 MR.It can be seen from Fig. 16, bottom right, however, that further penetration into the 10 MR channels by the ortho isomer leads to a very high increase in the interaction energy. Regarding diÜusion paths, we –rst note that in zeolites with parallel and perpendicular systems, visualisation of the channel structure is straightforward, but in a zeolite likeC. R.A. Catlow et al. 27 Fig. 16 Activation energies calculated for the diÜusion of ortho- and para-xylene through the channels of CIT-1 CIT-1 it is more difficult to understand the trajectories of the molecules in the channel system. In order to make the visualisation as clear as possible we have projected both the xylene trajectories and the channel axes on two-dimensional diagrams. The corresponding results for the three simulations are shown in Fig. 17, where only one of the three possible projections (xy, yz, xz) is given which is sufficient for our purpose. The three –gures correspond to the three diÜerent loadings given above. It can be seen [Fig. 17(a) and (b)] that ortho xylene diÜuses exclusively through the 12 MR channel system in both cases regardless of the loading simulated.Some penetration into the 10 MR channels can also be noted as was suggested by the calculated activation energies. The main feature of the diÜusion of ortho-xylene is that it tends to diÜuse in a restricted part of the zeolite with what we denote as extensive local motion. This behaviour is attributable to the high size of the ortho isomer with respect to the channel dimensions which precludes this molecule from travelling more freely through the structure.In the third simulation it can be seen [Fig. 17(c)] that para-xylene can diÜuse through both channel systems. The –gure also shows two additional interesting features : diÜusion through the 12 MR channel seems to be more favourable, as the higher population of this channel indicates ; and diÜusion through the 10 MR channel follows a straight line because the para isomer matches the channel size which precludes the para-xylene from rotating in this channel. Self-diÜusivities have been estimated from the history –les by calculating the mean square displacements of the centre of mass of the xylene molecules in the 100 ps simulation run.The results are illustrated and presented in Fig. 18 from which two main conclusions follow.The –rst concerns the in—uence of loading which can be estimated from the –rst two values corresponding to ortho-xylene. The results show that self-28 Computer modelling as a technique in solid state chemistry Fig. 17 Trajectory plots showing the diÜusion path of the centre of mass of the sorbate molecules in the CIT-1 macrocell: (a) 32 ortho-xylene, (b) 8 ortho-xylene, (c) 8 para-xyleneC.R. A. Catlow et al. 29 Fig. 18 Trajectories and diÜusion coefficients for the diÜerent simulated loadings and trajectory plot (100 ps) showing preferential diÜusion path through 12 MR channels. ortho-Xylene, 1.00 molecules uc~1, D\3.56]10~6 cm2 s~1; ortho-xylene, 0.25 molecules uc~1, D\7.79]10~6 cm2 s~1; para-xylene, 0.25 molecules uc~1, D\25.18]10~6 cm2 s~1.(a) View parallel to [001] (12 MR channels) ; (b) view parallel to [110] (10 MR channels). Trajectories of the centre of mass of the xylene molecule are shown in black. diÜusivities decrease with increasing loading, a feature common to many other zeolites due to the increasing importance of the repulsive guest»guest interactions.On the other hand, for the same loading we see a much higher self-diÜusivity for para-xylene due to its smaller activation energy and to the possibility of diÜusing also through the 10 MR channels. The present simulations have thus yielded semi-quantitative self-diÜusivities, and have yielded a good understanding of the microscopic features of the diÜusion process ; such results will be especially useful when combined with experiments in achieving a better understanding of diÜusion in zeolites. 7 Computer modelling and hydrothermal synthesis One of the most exciting recent developments in the –eld has been the exploitation of computational techniques in understanding the factors controlling the synthesis of microporous materials.A detailed discussion is given by Lewis et al.112 in this volume. Here we concentrate on one aspect of our work in this –eld, relating to the mechanisms of condensation of silica fragments»possibly the most fundamental process in hydrothermal synthesis.Condensation reactions of the type have been studied by Burggraf et al.,113 using semiempirical methods, and the activation energy is estimated experimentally to be 15 kcal mol~1.114 The reverse reaction has been studied by Lasaga and Gibbs115 (to explain the surface dissolution of quartz), who predicted an activation energy of 21.9 kcal mol~1, using ab initio calculations.The condensation reaction is discussed in detail by Iler116 and a comprehensive review is given by Brinker and Scherer.11730 Computer modelling as a technique in solid state chemistry In this work we studied the simplest condensation reaction of two species 2Si(OH)4 to form a dimer, using density functional theory (DFT) coupled with a continuum dielectric model (COSMO) to try to describe the electrostatic conditions found in real silica solutions. We used the Dmol code from MSI,118 with non-local (Becke»Lee» Yang»Parr)147,148 exchange and correlation energy and a DNP double numerical basis set. All atomic arrangements were –rst optimised in the gas-phase and subsequently recalculated with COSMO as a single energy point, without reoptimisation.Acid catalysis was considered throughout the work, corresponding to pH\3, the conditions mostly found in current experimental work. The protonated species are –rst investigated.The Si(OH)3(H2O)`, Si2 OH(OH)6 ` transformation of the protonated reactants into the protonated products is then analysed considering both an and a lateral attack mechanism. Assuming that the SN2 reaction occurs by the attack of a monomer on a second protonated monomer, the proton remaining after the water molecule has left should be attached to the oxygen that initiated the nucleophilic attack, i.e., the bridging oxygen.The reactions to protonate the reactant and the product can then be written as follows, with Si(OH)4 Si2OH(OH)6 ` the calculated energies (*E) given above the arrow: Si(OH)4]H3O` »»»»»»’ 7.4 kcal mol~1 Si(OH)3(H2O)`]H2O (III) Si2O(OH)6]H3O` »»»»»»’ 12.9 kcal mol~1 Si2OH(OH)6 `]H2O (IV) The calculated protonation energies agree with the basicity trends reviewed by Brinker and Scherer.117 The additional charge is more stabilised in the ion than in the H3O` protonated silicate clusters. The protonation energy is smaller for than for Si(OH)4 because the electron-withdrawing eÜect is smaller for OH than for OSi Si2OH(OH)6 ` groups, making the protonated monomer more stable.Assuming that these trends apply for larger clusters, in acid conditions small clusters and chain ends are more likely to be protonated and receive the nucleophilic attack of groups situated in the middle of the chains, forming branched chains.The reaction of the protonated monomer to give the protonated dimer is therefore an endothermic process, given by: Si(OH)4]Si(OH)3(H2O)` »»»»»»’ 5.1 kcal mol~1 Si2OH(OH)6 `]H2O (V) Two mechanisms for this reaction were studied considering both an and a lateral SN2 attack, represented in Fig. 19(a) and (b) respectively. The energy evolution along the reaction path shows –rst a small energy barrier SN2 of 2.5 kcal mol~1 before decreasing to a –ve-silicon intermediate, which is 3.4 kcal mol~1 more stable than the reactants. A second energy barrier of 11.3 kcal mol~1, the largest in the whole reaction process, occurs later, when the water molecule is leaving. In the lateral attack mechanism, after the –rst 2.5 kcal mol~1 energy barrier, a –vesilicon intermediate occurs, 1.9 kcal mol~1 less stable than the reactants. As in the gas phase, a pronounced peak occurs when the water molecule leaves, forming the largest energy barrier in the overall mechanism: 17.3 kcal mol~1, which is 6 kcal mol~1 larger than in the mechanism.SN2 The global energy evolution from the neutral reactants to the neutral products is depicted in Fig. 20. The activation energies calculated in this work for the (11.3 kcal SN2 mol~1) and lateral attack (17.3 kcal mol~1) mechanisms are in good agreement with the value estimated from experiment (15 kcal mol~1).Owing to the strategy utilised in this study, a gas-phase optimisation followed by a COSMO energy calculation, the diÜerence in energy between reactants and products is very small : [0.4 kcal mol~1. OnC. R. A. Catlow et al. 31 (a) and lateral attack (b) mechanisms for the condensation reaction (V) (see text). The Fig. 19 SN2 x axis shows the distance between the Si atoms of the dimer and the O of the (Aé ) Si2 OH(OH)6 ` departing The y axis shows the distance between the attacking O of the mol- H2O.(Aé ) Si(OH)4 ecule and the Si of the attacked monomer Si(OH)4H`. optimising both reactants and products in the COSMO environment, this diÜerence increases to [3.2 kcal mol~1. These results show that both the condensation and the hydrolysis reactions have a small activation energy.For both reactions, the mechanism seems to be a low SN2 energy route for converting the reactants into the products. However, as the results for the lateral attack show, several other mechanisms that are energetically and statistically less favourable should still be possible and may occur simultaneously in the solution.32 Computer modelling as a technique in solid state chemistry Fig. 20 Energy evolution (kcal mol~1) during the condensation reaction, for both Si(OH)4 SN2 and lateral attack mechanisms, from non-local DFT/COSMO calculations The DFT»COSMO implementation clearly oÜers a powerful tool to study the mechanisms occurring in this key condensation reaction. 8 Reactivity In this –nal section of recent applications we consider the use of quantum mechanical techniques in elucidating reaction mechanisms.Two topical examples are discussed. 8.1 The methanol to gasoline reaction The methanol to gasoline process119 in medium pore-sized aluminosilicate zeolites has been the subject of great theoretical120h133 and experimental119,134h146 interest due to its potential industrial importance but also, importantly, because of its status as a prototypical solid state acid catalysed reaction.Indeed, the methanol to gasoline Br‘nsted process is now used as a standard test for the presence of acidity. Despite the Br‘nsted obvious importance of this process we still have limited knowledge of the precise nature of the reaction, particularly with respect to the mechanism for formation of the –rst CwC bond.The past few years has seen a rapid increase in the use of DFT methods to study the electronic structure of chemically relevant systems, particularly since the introduction of accurate non-local corrections147h149 to the local exchange-correlation potentials150,151 that have been so successful in solid state physics. However, there are still questions concerning the validity of such methods in probing transition states of chemical reactions.This section therefore describes the study of a possible pathway for formation of the –rst CwC bond in the methanol to gasoline process with a number of quantum chemical methods in order to understand better the methanol to gasoline process and the limitations of the theoretical methods.Aluminosilicate zeolites consist of corner sharing tetrahedral units where T is TO4 either Al or Si. In the case of T\Al, a charge imbalance of [1 is generated which is compensated for by incorporation of a proton, H`, generating acidity. For Br‘nsted computational tractability only a few of these T atoms can be included in a quantum chemical calculation, a typical cluster model of the active site being a 3T fragment with dangling bonds terminated with hydrogen to preserve charge neutrality.However, a number of authors have used the smaller 1T model for a number of problems associatedC. R. A. Catlow et al. 33 with adsorption and reactivity at the acid site and found that results are sur- Br‘nsted prisingly similar to those obtained from a 3T model.127,152,153 One reason for this result could be that a typical process at an aluminosilicate acid site involves a push» Br‘nsted pull interaction with two of the oxygens attached to Al and although the 1T model is more basic at both oxygens than the 3T model, the diÜerence in the basic properties of the two oxygens in the 1T and 3T models is expected, due to symmetry, to be similar.However, it must be remembered that the electrostatic potential above the 1T model is more negative than that of the 3T model and the latter would therefore stabilise the localisation of negative charge more than the former, and vice versa. We therefore employ a hydrogen terminated 1T model, or ZOH, and (HO)3Al(OH2) its methylated analogue, or ZOMe, to study the formation of an ethyl- (HO)3AlO(Me)H ated surface, or ZOEt, and from that ethene, see Fig. 21. Formation of (HO)3AlO(Et)H the surface methoxy, ZOMe, has been the subject of previous studies123,127,129,130 and will not be considered further except to note that it can be formed at the site Br‘nsted by two methanol molecules with an activation barrier of about 130 kJ mol~1. We have studied this reaction at the HF/6-31G**//HF/6-31G**, the MP2/6»31G**//HF/6- 31G**, the BLYP/DZVP//LDA/DZVP and the BLYP/DZVP//BLYP/DZVP levels of theory (notation indicates the level of theory used for –nal energy calculations//levels of theory used for geometry optimisations).The Hartree»Fock and MP2 calculations were performed using the GAUSSIAN94 code.154 The local density approximation (LDA)150,151 and the gradient corrected Becke147»Lee»Yang»Parr148 (BLYP) DFT calculations were performed with the Cray Research code DGauss 3.0155 using the HF/6- 31G** optimised structures as starting points. Table 8 summarises the results.All structures were similar to the HF/6-31G** structures shown in Fig. 21. Fig. 21 HF/6-31G** reaction path for formation of ethene from a surface methoxy (ZOMe) and adsorbed methanol (MeOH) via a carbene-like transition state (ts-I) and surface ethoxy intermediates (int-I and int-II).The acid site was modelled by a 1T, fragment. Br‘nsted (OH)3AlO(H)H The energy changes are given in Table 8.34 Computer modelling as a technique in solid state chemistry Table 8 Energy changes (kJ mol~1) to Fig. 21. methoda *E(1)b *E(2) *E(3) *E(4)b *E(5) HF/6-31G**//HF/6-31G** 514 [573 39 232 [172 MP2/6-31G**//HF/6-31G** 465 [539 53 194 [126 BLYP/DZVP//LDA/DZVP 386 [414 86 136 [78 BLYP/DZVP//BLYP/DZVP 382 [419 37 140 [88 a Notation indicates method/basis used for calculating the –nal energy//method/ basis used for optimising the geometry.b Activation energy barrier. The reaction path shown in Fig. 21 describes the deprotonation of a surface methoxy group leading to CwC bond formation (in ZOEt) via a transition state with a surface stabilised carbene-like species (ts-I).Fracture of the CwH bond occurs early in the reaction coordinate to form water with the OH of methanol and this bond fracture is probably the main contributor to the activation barrier. Deprotonation of this ZOEt species results in the formation of adsorbed ethene through the transition state ts-II.It is immediately clear that the diÜerent quantum chemical approximations lead to quite diÜerent estimates of the energy changes, especially for the –rst stage of the reaction. It is expected that the HF//HF results should overestimate the activation barriers as the adiabatic transition states are generally formed from an avoided crossing of the ground and one or more excited states156 and therefore the single determinantal HF result for the transition state is expected to be a poorer approximation to the solution of the Schroé dinger equation than are the HF results for the minimum energy structures.Similarly, the MP2//HF activation barrier, although better than the HF//HF result, is still probably an overestimate as the MP2 correction will not recover all of the correlation energy.The BLYP//LDA and BLYP//BLYP results are somewhat surprising as the geometries and the energy changes are remarkably similar (all structures are within 0.2 and 10° of each other) although the calculated activation barriers are prob- ” ably still underestimates.157 This is in contrast to a number of reports suggesting that LDA is particularly poor at reproducing hydrogen bonds (see, for example, ref. 130). Uncertainties in the model and method aside, it is clear that deprotonation of surface methoxy species to form ZOEt will be unlikely, although, as discussed above, use of a larger model will result in a less negative electrostatic potential above the cluster and therefore in a stabilisation of the carbene-like, electron-rich transition state, ts-I.As a comparison with other work, Blaszkowski and van Santen132 have studied a number of reaction pathways with a 1T model at the BP (Becke147»Perdew149)//LDA level of theory in which CwC bond formation occurs in the adsorbate, not the surface carbon species. The resulting calculated activation barriers are substantially lower than those presented here.Calculations on larger clusters will be reported in future work. 8.2 Partial oxidation on oxide surfaces Facilitating CwH bond breakage is one of the typical functions of catalysts. The oxidative coupling of methane (OCM) represents an important process, the initial step of which requires just that : methyl radicals are formed from methane by H-abstraction.One of the factors determining the usefulness of a catalyst for such a process will be the amount by which the energy barrier for the initial CwH bond cleavage is reduced in its presence. A variety of oxides of diÜerent complexity are known to act as OCM catalysts ;158 in some cases their catalytic activity can be promoted by dopants.It has been suggested that these promoters play a role in stabilising electronic defect sites on the catalysts surface.159 For example, localisation of electron hole states preferentially takes place atC. R. A. Catlow et al. 35 anion surface sites adjacent to cation substitutions, where the ionic charge of the dopant is lower than that of the original cation ; the resulting O~ species (known to be present from EPR measurements160) may very well function as active sites in hydrogen abstraction reactions.It has long been proposed that [Li]0 centres constitute such a catalytically active site on the surface of Li doped magnesium oxide surfaces. This would explain why doping is found to have an enhancing eÜect on this very extensively studied OCM catalyst.161 The [Li]0 defect site, consisting of a neighbouring pair of a Li` (substituting for Mg2`) and an electron hole localised on an oxygen ion, is thought to interact with the hydrocarbon molecule in the following process : CH4][Li`O~]]~CH3][Li`(OH)~] (VI) leading to the gas phase radical species observed experimentally.162 ~CH3 Despite continued interest in this system (see, for example, ref. 158 and 161) little is known about reaction (VI) in detail. A number of theoretical studies have been devoted to this question (and related ones) during the last ten years163h176 and even further back.177h179 Most of them are performed at the HF level164,172,177h179 or below.163,165,170,174 On post-HF level studies involving MP2,170,172,174,180 CASSCF,167h169 DFT175 and hybrid176 approaches seem to indicate that electron correlation may well play a signi–cant role in the proper description of the system. Although the eÜect of Li doping is mentioned in most of these studies, only a few incorporate it explicitly in the models used.168,169,174,176,180 In some, relaxation of the substrate is taken into account178h180 but in other cases it is neglected.Finally some of the references deal only with dissociation, not with methane.H2 Taking all this into account it seems desirable to obtain a theoretical description of this system using (a) an explicit model for the [Li]0»methane interaction on a non-rigid MgO and employing (b) a theoretical approach that allows us to estimate the importance of electron correlation in this reaction.It is not easy, however, to satisfy both goals simultaneously at reasonable computational cost. A cluster model of considerable size is necessary to represent substrate relaxation taking place in the course of reaction (VI). On the other hand, to assess whether ab initio calculations at the HF level are sufficient to describe reaction (VI) it is necessary to test ways of including electron correlation to some degree.Using MP2 or MP4 corrections to the HF approach, one is limited to cluster sizes that do not comply with (a) above. Therefore, it is important to validate less expensive alternatives such as, for example, DFT. This has been done, using a very simple model system in which the [Li]0 centre is represented by a minimal cluster Li`O~ embedded in a point charge array.181 Comparing the energy pro–le along the reaction coordinate (Fig. 22) one observes the eÜect of electron correlation : from HF calculations a symmetric energy barrier is obtained; the height of this barrier is substantially reduced as MP2 corrections are added (even more for MP4). DFT (using BLYP gradient corrections147,148) on the other hand reduces the activation barrier to a mere shoulder in the corresponding curve.For future studies employing larger clusters we therefore conclude that while results obtained from HF calculations may be expected to overestimate the barrier, those coming from DFT calculations tend give too low a value and should be used with caution. The next step in reaching a realistic description of reaction (VI) requires the construction of suitable cluster models.To represent ionic relaxation adequately, more than only the nearest neighbours of the active site will have to be included (here it has to be taken into account that the outermost atoms of the cluster must be kept –xed in relation to the surrounding point charges). Enlarging the cluster, however, not only con—icts with computational economy but also leads to SCF convergence problems, especially in the case of DFT.A cluster of the following composition has been used successfully : –rst layer (12 Mg –xed); second layer (8 Mg –xed); third layer (4 Mg16O9 Mg8LiO4 Mg4O36 Computer modelling as a technique in solid state chemistry Fig. 22 Energy pro–le of reaction (VI) as calculated from HF (»»), MP2 (- - -), MP4 (- … - …-) and DFT»BLYP (…………).The sum of the energies of the products in reaction (VI) has been used as an energy oÜset for each of these curves. Mg –xed) and fourth layer Mg (–xed). Preliminary calculations of the [Li]0 centre yield a relaxation energy of 1.9 eV with ionic displacements of up to 0.3 (compared to the ” bulk positions). At the present stage in this study it is clear that both aspects discussed above play a crucial role : (a) the quality of the model (including geometric relaxation) and (b) the quality of the method (beyond HF).Consequently we shall combine the two approaches described here in our future work, leading to a more accurate determination of CwH bond activation by Li/MgO. Further details of the results described above are given in ref. 181. 9 Conclusions The applications described in this article illustrate, we hope, the range and diversity of the growing –eld of computational solid state chemistry. We can expect a continuing increase in the complexity of the systems studied and the accuracy of the predictions made, enhancing the value of computational tools as techniques in solid state chemistry. We thank EPSRC for supporting this work. We also gratefully acknowledge –nancial support from DSM, ICI, and MSI.We are grateful to W. C. Mackrodt, A. M. Stoneham, A. K. Cheetham and Sir John Meurig Thomas for valuable discussions on many aspects of this work. References 1 C. R. A. Catlow, R. G. Bell and J. D. Gale, J. Mater. Chem., 1994, 4, 781. 2 M. P. Allen and D. J. Tildesley, in Computer Simulation of L iquids, Oxford University Press, 1987. 3 C. R. A. Catlow and G. D. Price, Nature (L ondon), 1990, 347, 243. 4 N. L. Allan and W. C. Mackrodt, Adv. Solid State Chem., 1993, 3, 221.C. R. A. Catlow et al. 37 5 Computer Modelling in Inorganic Crystallography, ed. C. R. A. Catlow, Academic Press, London, 1997. 6 C. R. A. Catlow, J. M. Thomas, C. M. Freeman, P.A. Wright and R. G. Bell, Proc. R. Soc. L ondon A, 1993, 442, 85. 7 J. H. Harding, Rep. Prog. Phys., 1990, 53, 1403. 8 N. F. Mott and M. J. Littleton, T rans. Faraday Soc., 1938, 34, 389. 9 S. C. Parker and G. D. Price, Adv. Solid State Chem., 1996, 1, 295. 10 M. J. Gillan, Physica B, 1985, 131, 157. 11 M. T. Dove and G. S. Pawley, J. Phys. C: Solid State Phys., 1984, 17, 6581. 12 R. W. Impey, M. Sprik and M. L. Klein, J. Chem. Phys., 1985, 83, 3638. 13 S. Yashonath, J. M. Thomas, A. K. Nowak and A. K. Cheetham, Nature (L ondon), 1988, 331, 601. 14 C. M. Freeman, C. R. A. Catlow, J. M. Thomas and S. Chem. Phys. L ett., 1991, 186, 137. Br‘de, 15 A. Shubin, C. R. A. Catlow, J. M. Thomas and K. I. Zamaraev, Proc. R. Soc. L ondon A, 1994, 446, 411. 16 G. E. Murch, Philos.Mag. A, 1982, 46, 151. 17 G. E. Murch, Philos. Mag. A, 1982, 46, 575. 18 A. D. Murray, G. E. Murch and C. R. A. Catlow, Solid State Ionics, 1985, 18/19, 196. 19 S. M. Auerbach, L. M. Bull, N. J. Henson, H. L. Metu and A. K. Cheetham, J. Phys. Chem., 1996, 100, 5923. 20 Computer simulations of solids, L ecture Notes in Physics, ed. C. R. A. Catlow and W. C. Mackrodt, Springer-Verlag, 1982, vol. 166. 21 Special Issue : Interatomic Potentials, Philos. Mag. B, 1996, 7. 22 B. G. Dick and A. W. Overhauser, Phys. Rev. B, 1958, 112, 90. 23 P. A. Madden and M. Wilson, Chem. Soc. Rev., 1996, 25, 339. 24 M. Wilson, P. A. Madden and B. J. Costa-Cabral, J. Phys. Chem., 1996, 100, 1227. 25 N. M. Harrison and M. Leslie, Mol. Simul., 1992, 9, 171. 26 J. D. Gale, C. R. A.Catlow and W. C. Mackrodt, Modell. Simul. Mater. Sci. Eng., 1992, 1, 73. 27 A. Takada, C. R. A. Catlow and G. D. Price, J. Phys. C: Solid State Phys., 1995, 7, 8659. 28 R. Car and M. Parrinello, Phys. Rev. L ett., 1985, 55, 2471. 29 R. Dovesi, C. Roetti, C. Freyria-Fava, E. Apra, V. R. Saunders and N. M. Harrison, Philos. T rans. A, 1992, 341, 203. 30 J. C. White and A. C. Hess, J.Phys. Chem., 1993, 97, 8703. 31 M. Sierka and J. Sauer, Faraday Discuss., 1997, 106, 41. 32 P. Sherwood, A. H. deVries, S. J. Collins, S. P. Greatbanks, N. A. Burton, M. A. Vincent and I. H. Hillier, Faraday Discuss., 1997, 106, 79. 33 A. M. Gorman, C. M. Freeman, C. M. Koé lmel and J. M. Newsam, Faraday Discuss., 1997, 106, 489. 34 C. R. A. Catlow, in Handbook of Catalysis, ed.G. Ertl, H. Knoé zinger and J. Weitcamp, Wiley-VCH, Weinheim, 1997. 35 W. R. Busing, WMIN, A computer program to model molecules and crystals in terms of potential energy functions, ORNL-5747, Oak Ridge National Laboratory, Oak Ridge, 1981. 36 C. R. A. Catlow, A. N. Cormack and F. Theobald, Acta Crystallogr., Sect. B, 1984, 40, 195. 37 M. Leslie, Daresbury Laboratory, Warrington, Cheshire, UK WA4 4AD. 38 J. D. Gale, J. Chem. Soc., Faraday T rans., 1997, 93, 629. 39 R. G. Bell, R. A. Jackson and C. R. A. Catlow, J. Chem. Soc., Chem. Commun, 1990, 782. 40 M. J. Henson, A. K. Cheetham and J. D. Gale, Chem. Mater., 1994, 6, 147. 41 M. D. Shannon, J. L. Cusci, P. A. Cox and S. J. Andrews, Nature (L ondon), 1991, 353, 417. 42 T. S. Bush, C. R. A. Catlow and P. D. Battle, J.Mater. Chem., 1995, 5, 1269. 43 C. T. Kresge, M. E. Leonowicz, W. J. Roth, J. C. Vartuli and J. S. Beck, Nature (L ondon), 1992, 359, 710. 44 J. S. Beck, J. C. Vartuli, W. J. Roth, M. E. Leonowicz, C. T. Kresge, K. D. Schmitt, C. T-W. Chu, D. H. Olsen, E. W. Sheppard, S. B. McCullen, J. B. Higgins and J. L. Schlenker, J. Am. Chem. Soc., 1992, 114, 10834. 45 P. T. Tanev, M. Chibwe and T.J. Pinnavaia, Nature (L ondon), 1994, 368, 321. 46 F. Rey, G. Sankar, T. Maschmeyer, J. M. Thomas, R. G. Bell and G. N. Greaves, T op. Catal., 1996, 3, 121. 47 A. Jentys and C. R. A. Catlow, Catal. L ett., 1993, 22, 251. 48 J. V. Smith and W. J. Dytrych, Nature (L ondon), 1984, 309, 607. 49 A. Navrotsky, I. Petrovic, Y. T. Hu, C. Y. Chen and M. E. Davis, Microporous Mater., 1995, 4, 95. 50 P. A. Cox, T he electronic structure and chemistry of solids, Oxford University Press, 1987. 51 J. K. Burdett, Chemical bonding in solids, Oxford University Press, New York, 1995. 52 C. Pisani, R. Dovesi and C. Roetti, Hartree»Fock Ab Initio T reatment of Crystalline Systems (L ecture Notes in Chemistry, V ol. 48), Springer, Heidelberg, 1988; 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, Italy, 1996. 53 M. Methfessel, Phys. Rev. B, 1988, 38, 1537; M. Methfessel, C. O. Rodriguez and O. K. Anderson, Phys. Rev. B, 1989, 40, 2009.38 Computer modelling as a technique in solid state chemistry 54 M. G. Stachiotti, F. Cora` , C. R. A. Catlow and C.O. Rodriguez, Phys. Rev. B, 1997, 55, 7508; F. Cora` , M. G. Stachiotti, C. R. A. Catlow and C. O. Rodriguez, J. Phys. Chem. B, 1997, 101, 3945. 55 R. A. Wheeler, M-H. Whangbo, T. Hughbanks, R. HoÜman, J. K. Burdett and T. A. Albright, J. Am. Chem. Soc., 1986, 108, 2222. 56 J. K. Burdett and S. A. Gramsch, Inorg. Chem., 1994, 33, 4309. 57 F. Cora` , A. Patel, N. M. Harrison, R. Dovesi and C.R. A. Catlow, J. Am. Chem. Soc., 1996, 118, 12174. 58 F. Cora` , A. Patel, N. M. Harrison, C. Roetti and C. R. A. Catlow, J. Mater. Chem., 1997, 7, 959. 59 D. H. Gay and A. L. Rohl, J. Chem. Soc., Faraday T rans., 1995, 91, 925. 60 P. W. Tasker, Philos. Mag. A, 1979, 39, 119; AERE Harwell Report, 1978, R9130. 61 J. G. Fripiat, A. A. Lucas, J. M. Andreç and E. G. Derouane, Chem.Phys., 1977, 21, 101. 62 W. C. Mackrodt, R. J. Davey, S. N. Black and R. Docherty, J. Cryst. Growth, 1987, 80, 441. 63 W. C. Mackrodt, Philos. T rans. R. Soc. L ondon, 1992, 341, 301. 64 P. Guenard, G. Renaud, A. Barbier and M. Gautier-Soyer, Mat. Res. Soc. Symp. Proc., 1996, 437, 15. 65 C. R. A. Catlow, J. Phys. Chem. Solids., 1977, 38, 1131, and references therein. 66 M.A. Nygren, D. H. Gay and C. R. A. Catlow, Surf. Sci., 1997, 380, 113. 67 N. H. de Leeuw, G. W. Watson and S. C. Parker, J. Phys. Chem., 1995, 99, 17219 and references therein. 68 J. P. Perdew and Y. Wang, Phys. Rev. B, 1986, 33, 8800; J. P. Perdew, Phys. Rev. B., 1986, 33, 8832. The code used for the calculations was deMon-KS3p1 written by D. R. Salahub, R. Fournier, P. Mlynarski, I.Papai, A. St-Amant and J. Ushio, in Density Functional Methods in Chemistry, ed. Labanowski and J. Andzelm, Springer, New York, 1991, p. 77; A. St-Amant, PhD Thesis, Universiteç de Montreç al, 1992. The present version of the code has been substantially modi–ed by L. G. M. Pettersson. 69 F. Bertaut, Compt. Rend., 1958, 246, 3447. 70 F. H. Jones, K. Rawlings, J. S. Foord, R.G. Egdell, J. B. Pethica, B. M. R. Wanklyn, S. C. Parker and P. M. Oliver, Surf. Sci., 1996, 359, 107. 71 R. Dovesi, M. Causa, R. Orlando, C. Roetti and V. R. Saunders, J. Chem. Phys., 1990, 92, 7402. 72 R. Docherty, K. J. Roberts, V. R. Saunders, S. Black and R. J. Davey, Faraday Discuss., 1993, 95, 11. 73 G. Clydesdale, K. J. Roberts and R. Docherty, J. Cryst. Growth, 1996, 166, 78. 74 A. R. George, K. D. M. Harris, A. L. Rohl and D. H. Gay, J. Mater. Chem., 1995, 5, 133. 75 P. P. Ewald, Ann. Phys., 1921, 64, 253. 76 Z. Berovitch-Yellin, J. van Mil, L. Addadi, M. Idelson, M. Lahav and L. Leiserowitz, J. Am. Chem. Soc., 1985, 107, 3111. 77 Z. Berovitch-Yellin, J. Am. Chem. Soc., 1985, 107, 8239. 78 A. L. Rohl and D. H. Gay, J. Cryst. Growth, 1996, 166, 84. 79 N.L. Allan, A. L. Rohl, D. H. Gay, C. R. A. Catlow, R. J. Davey and W. C. Mackrodt, Faraday Discuss., 1993, 95, 273. 80 L. Verlet, Phys. Rev., 1967, 159, 98. 81 W. C. Swope, H. C. Anderson, P. H. Berens and K. R. Wilson, J. Phys. Chem., 1982, 76, 637. 82 C. W. Gear, in T he numerical integration of ordinary diÜerential equations of various orders, report ANL 7126, 1966. 83 C. W.Gear, in Numerical initial value problems in ordinary diÜerential equations, Prentice-Hall, Englewood CliÜs, NJ, 1971. 84 E. S. Boek, W. J. Brierls, J. van Eerden and D. Feil, J. Chem. Phys., 1992, 96, 7010. 85 A. T. Hagler, S. Lifson and P. Dauber, J. Am. Chem. Soc., 1979, 101, 5122. 86 Molecular Simulations, 9685 Scranton Road, San Diego, USA. 87 M. Hatada, J. Janacorick, B. Graves and S-H.Kim. J. Am. Chem. Soc., 1985, 107, 4279. 88 J. D. H. Donnay, Am. Mineral., 1937, 22, 446. 89 H. J. C. Berendsen, J. P. M. Postma, W. F. van Guntersen and J. Hermans, in Intermolecular Forces, ed. B. Pullman, Reidel, Dordrecht, 1981. 90 K. V. Wright and C. R. A. Catlow, Phys. Chem. Miner., 1994, 20, 515. 91 K. V. Wright and C. R. A. Catlow, Phys. Chem. Miner., 1996, 23, 38. 92 D. R. Bell and G. R. Rossman, Science, 1992, 255, 1391. 93 Q. Bai and D. L. Kohlstedt, Phys. Chem. Miner., 1993, 19, 460. 94 C. R. A. Catlow, P. S. Baram, S. C. Parker, J. Purton and K. V. Wright, Proc. on the 7th Europhysical Conference on Defects in Insulating Materials, Eurodim 94, Lyon (1994), 1995, 134, 57. 95 C. R. A. Catlow, P. S. Baram, S. C. Parker, J. Purton and K. V. Wright, Philos.T rans. R. Soc. L ondon A, 1995, 350, 265. 96 J. Rabo, in Zeolite Chemistry and Catalysis, ACS Monograph Series, vol. 171, 1976. 97 J. Scherzer, Catal. Rev. Sci. Eng., 1989, 31, 215. 98 A. Corma, Chem. Rev., 1995, 95, 559. 99 A. Navrotsky, J. Petrovic, Y. Hu, C-Y. Chen and M. E. Davis, Microporous Mater., 1995, 4, 95. 100 S. M. Csicsery, Zeolites, 1983, 4, 202.C.R. A. Catlow et al. 39 101 D. S. Santilli and S. J. Zones, Catal. L ett., 1990, 7, 383. 102 D. S. Santilli, T. V. Harris and S. I. Zones, Microporous Mater., 1993, 1, 329. 103 R. M. Barrer, in Zeolites and Clay Minerals, Academic Press, London, 1978. 104 J. Kaé rger and J. Caro, J. Chem. Soc., Faraday T rans., 1977, 1, 1363. 105 H. Jobic, M. Bee, J. Caro, M. Bulow and J. Kaé rger, J.Chem. Soc., Faraday T rans., 1989, 85, 4201. 106 (a) J. B. Nicholas, F. R. Trouw, J. E. Mertz, L. E. Iton and A. J. Hop–nger, J. Phys. Chem., 1993, 97, 4149; (b) H. Klein, H. Fuess and G. Schrimpf, J. Phys. Chem., 1996, 100, 11101; (c) S. M. Auerbach, N. J. Henson, A. K. Cheetham and H. I. Metiu, J. Phys. Chem., 1995, 99, 10600. 107 R. F. Lobo and M. E. Davis, J. Am. Chem.Soc., 1995, 117, 3766. 108 L. B. Young, S. A. Butter and W. W. Kaeding, J. Catal., 1982, 76, 418. 109 G. Sastre, N. Raj, C. R. A. Catlow, R. Roque-Malherbe and A. Corma, J. Phys. Chem., in press. 110 T. R. Forrester and W. Smith, in DLñPOLYñ2.0 User Manual, CLRC Daresbury Laboratory, 1995. 111 Catalysis6.0 Software Package, MSI San Diego, 1995. 112 D. W. Lewis, C. R. A. Catlow and J.M. Thomas, Faraday Discuss., 1997, 106, 451. 113 L. W. Burggraf, L. P. Davis and M. S. Gordon, in Ultra-structure Processing of Adanced Materials, ed. D. R. Uhlmann and D. R. Ulrich, John Wiley and Sons, New York, 1992. 114 S. H. Garofalini and G. Martin, J. Phys. C: Solid State Phys., 1994, 98, 1311. 115 A. C. Lasaga and G. V. Gibbs, Am. J. Sci., 1990, 290, 263. 116 R. K. Iler, in T he Chemistry of Silica, John Wiley and Sons, 1979. 117 C. J. Brinker and G. W. Scherer, Sol»Gel Science : T he Physics and Chemistry of Sol»Gel Processing, Academic Press, 1989. 118 DMOL manual, Molecular Simulations, 1996. 119 S. L. Meisel, J. P. McCullogh, C. H. Lechthaler and P. B. Weisz, CHEMT ECH 1976, 6, 86. 120 F. Haase and J. Sauer, J. Am. Chem. Soc., 1995, 117, 3780. 121 S.Greatbanks, I. H. Hillier, N. A. Burton and P. Sherwood, J. Chem. Phys., 1996, 105, 3770. 122 J. Limtrakul, Chem. Phys., 1995, 193, 79. 123 C. M. Zicovich-Wilson, P. Viruela and A. Corma, J. Phys. Chem., 1995, 99, 13224. 124 E. Nusterer, P. E. Bloé ch and K. Schwartz, Angew. Chem., Int. Ed. Engl., 1996, 35, 175. 125 E. Nusterer, P. E. Bloé chl and K. Schwartz, Chem. Phys. L ett., 1996, 253, 448. 126 R. Shah, M. C. Payne, M-H. Lee and J. D. Gale, Science, 1996, 271, 1395. 127 P. E. Sinclair and C. R. A. Catlow, J. Chem. Soc., Faraday T rans., 1996, 92, 2099. 128 P. E. Sinclair and C. R. A. Catlow, J. Chem. Soc., Faraday T rans., 1997, 93, 333. 129 P. E. Sinclair and C. R. A. Catlow, J. Phys. Chem., 1997, 101, 295. 130 S. R. Blaszkowski and R. A. van Santen, J. Phys. Chem., 1995, 99, 11728. 131 S. R. Blaszkowski and R. A. van Santen, J. Am. Chem. Soc., 1996, 118, 5152. 132 S. R. Blaszkowski, PhD Thesis, Technical University of Eindhoven, 1997 (supervised by Prof. R. van Santen). 133 R. Shah, J. D. Gale and M. C. Payne, J. Phys. Chem., submitted. 134 T. R. Forester and R. F. Howe, J. Am. Chem. Soc., 1987, 109, 5076. 135 U. Messow, K. Quitzsch and H. Herden, Zeolites, 1984, 4, 255. 136 C. D. Chang and A. J. Silvestri, J. Catal., 1977, 47, 249. 137 Y. Ono and T. Mori, J. Chem. Soc., Faraday T rans., 1981, 77, 2209. 138 J. P. van den Berg, J. P. Wolthuizen and J. H. C. van Hoof, Proc. 5th Int. Zeolite Conf., ed. L. V. C. Rees, Heydon, London, 1981, 649. 139 G. A. Olah, Pure Appl. Chem., 1981, 53, 201. 140 J. K. A. Clarke, R. Darcy, B. F. Hegarty, E. OœDonoghue, V. Ebrahimi and J. J. Rooney, J. Chem. Soc., Chem. Commun., 1986, 425. 141 G. J. Hutchings and R. Hunter, Catal. T oday, 1990, 6, 279 (review article»see also references therein). 142 J. E. Jackson and F. M. Bertsch, J. Am. Chem. Soc., 1990, 112, 9085. 143 J. Bandiera and C. Naccache, Appl. Catal., 1991, 69, 139. 144 E. J. Munson, A. A. Kheir and J. F. Haw, J. Phys. Chem., 1993, 97, 7321. 145 E. J. Munson, A. A. Kheir, N. D. Lazo and J. F. Haw, J. Phys. Chem., 1992, 96, 7740. 146 E. J. Munson and J. F. Haw, J. Am. Chem. Soc., 1991, 113, 6303. 147 A. D. Becke, Phys. Rev. A, 1988, 38, 3098. 148 C. Lee, W. Yang and R. G. Parr, Phys. Rev. B, 1988, 37, 785. 149 J. P. Perdew, Phys. Rev. B, 1986, 33, 8822. 150 S. H. Vosko, L. Wilk and M. Nusair, Can. J. Phys., 1980, 58, 1200. 151 P. A. M. Dirac, Proc. Cambridge Philos. Soc., 1930, 26, 376. 152 E. M. Evleth, E. Kassab and L. R. Sierra, J. Phys. Chem., 1994, 98, 1421. 153 E. M. Evleth, E. Kassab, H. Jessri, M. Allavena, L. Montero and L. R. Sierra, J. Phys. Chem., 1996, 100, 11368. 154 GAUSSIAN94, rev. D, Gaussian, Pittsburgh, PA, 1994. M. J. Frisch, G. W. Trucks, H. B. Schlegel, P. M. W. Gill, M. A. Cheeseman, J. R. Robb, T. A. Keith, G. A. Petersson, J. A. Montgomery, K. Raghavachari, M. A. Al-Laham, V. G. Zakrzewski, J. V. Ortiz, J. B. Foresman, J. Cioslowski, B. B.40 Computer modelling as a technique in solid state chemistry Stefanov, N. Nanayakkara, M. Challacombe, C. Y. Peng, P. Y. Ayala, W. Chen, M. W. Wong, J. L. Andreç s, E. S. Replogle, R. Gomperts, R. L. Martin, D. J. Fox, J. S. Binkley, D. J. Defrees, J. Baker, J. J. P. Stewart, M. Head-Gordon, C. Gonzalez, J. A. Pople. 155 J. Andzelm and E. Wimmer, J. Phys. Chem., 1992, 96, 1280. DGauss 3.0 is a density functional theory electronic structure code from Cray Research. 156 S. S. Shaik, H. B. Schlegel and S. Wolfe, in T heoretical Aspects of Physical Organic Chemistry: T he Mechansim, J. Wiley and Sons, New York, 1992. SN2 157 L. Fan and T. Ziegler, J. Am. Chem. Soc., 1992, 114, 10890. 158 G. J. Hutchings, M. S. Scurrell and J. R. Woodhouse, Chem. Soc. Rev., 1989, 18, 251. 159 T. Ito, J-X. Wang, C-H. Lin and J. H. Lunsford, J. Am. Chem. Soc., 1985, 107, 5062. 160 T. Ito and J. H. Lunsford, Nature (L ondon), 1985, 314, 721. 161 J. H. Lunsford, Angew. Chem., Int. Ed. Engl., 1995, 34, 970. 162 D. J. Driscoll, W. Martir, J-X. Wang and J. H. Lunsford, J. Am. Chem. Soc., 1985, 107, 58. 163 S. P. Mehandru, A. B. Anderson and J. F. Brazdil, J. Am. Chem. Soc., 1988, 110, 1715. 164 C. M. Zicovich-Wilson, R. Gonzaç lez-Luque and P. M. Viruela-Martiç n, J. Mol. Struct. (T HEOCHEM), 1990, 208, 153. 165 N. U. Zhanpeisov, A. G. Pelœenschchikov and G. M. Zhidomirov, Kinet, Catal. (Engl. T rans.), 1990, 31, 491. 166 H. Kobayashi, M. Yamaguchi and T. Ito, J. Phys. Chem. 1990, 94, 7206. 167 K. J. J. Chem. Phys., 1991, 95, 4626. B‘rve, 168 K. J. and L. G. M. Petersson, J. Phys. Chem., 1991, 95, 3214. B‘rve 169 K. J. Boé rve and L. G. M. Pettersson, J. Phys. Chem., 1991, 95, 7401. 170 P. M. Viruela-Martiç n, R. Viruela-Martiç n, C. M. Zicovich-Wilson and T. Tomaç s-Vert, J. Mol. Catal., 1991, 64, 191. 171 T. Ito, T. Tashiro, M. Kawasaki, T. Watanabe, K. Toi and H. Kobayashi, J. Phys. Chem., 1991, 95, 4476. 172 K. Sawabe, N. Koga, K. Morokuma and Y. Iwasawa, J. Chem. Phys., 1993, 97, 6871. 173 J. L. Anchell, K. Morokuma and A. C. Hess, J. Chem. Phys., 1993, 99, 6004. 174 M-a. D. Stiakaki, A. C. Tsipis, C. A. Tsipis and C. E. Xanthopoulos, J. Mol. Catal., 1993, 82, 425. 175 H. Kobayashi, D. R. Salahub and T. Ito, J. Phys. Chem., 1994, 98, 5487. 176 R. Orlando, F. Cora` , R. Millini, G. Perego and R. Dovesi, J. Chem. Phys., 1996, 105, 8937. 177 E. G. Derouane, J. G. Fripiat and J. M. Andreç , Chem. Phys. L ett., 1974, 28, 445. 178 E. A. Colbourn and W. C. Mackrodt, Surf. Sci., 1982, 117, 571. 179 S. A. Pope, M. F. Guest, I. H. Hillier, E. A. Colbourn, W. C. Mackrodt and J. Kendrick, Phys. Rev. B, 1983, 28, 219. 180 M. A. Johnson, E. V. Stefanovich and T. N. Truong, J. Phys. Chem. B, 1997, 101, 3196. 181 L. Ackermann, J. D. Gale and C. R. A. Catlow, J. Phys. Chem., in press. Paper 7/01964E; Received 20th March, 1997