C Springer-Verlag 1997 Phys Chem Minerals (1997) 25: 70–78 OR IG I NA L P A P E R G.W. Watson ? P.M. Oliver ? S.C. Parker Computer simulation of the structure and stability of forsterite surfaces Received: 28 October 1996 y Revised, accepted: 7 April 1997 Abstract The aim of this paper is to demonstrate that atomistic simulations can be used to evaluate the structure of mineral surfaces and to provide reliable data for forsterite surfaces up to a plane index of 2 using the code METADISE. The methods used to calculate the surface structure and energy which have more commonly been used to study ceramics are briefly explained as is a comparison with experimental data, most notable the crystal morphology. The predicted morphologies show that all the methods (Donnay-Harker, Attachment energies and equilibrium) show most of the surfaces that are expressed in observed crystals. The equilibrium morphology calculated from the relaxed surface energies is the only method which expresses the {201} surfaces and the {101} surfaces, which appear only upon relaxation. The more stable surfaces are shown to be those which have the highest surface density and more closely resemble close packed structures with highly coordinated surface ions and silicon as far from the surface as possible. The most stable surfaces the {100} which has alternating layers of MgO and SiO2 terminating with an MgO layer. The structure is similar to the MgO {100} surfaces and has a similar energy (1.28 Jm22 compared to 1.20). The second most stable are the {201} which have a stepped surface topology, but is also compact with a relaxed surface energy of 1.56 Jm22. The results indicate that atomistic simulation is well suited to the prediction of surface structure and morphology although care must be taken in choosing potentials which model the structure and elastic properties accurately. G.W. Watson1 (✉) ? P.M. Oliver ? S.C. Parker Computational Solid State Chemistry Group, School of Chemistry, University of Bath, Claverton Down, Bath, BA2 7AY. U.K. Fax: 144(0)1225 826 231 email: g.watson @liv.ac.uk Present address: Leverhulme Centre for Innovative Catalysis, Department of Chemistry, University of Liverpool, Oxford Road, Liverpool, L69 3BX, U.K. Introduction Atomistic simulation of the surfaces of polar solids was pioneered by the work of Tasker (1978, 1979a) and Mackrodt and Stewart (1977). Initially, these simulations concentrated on simple ceramics such as the cubic halides (Tasker 1979b) and binary oxides, MgO (Colbourn et al. 1983) and NiO (Tasker and Duffy 1984). Recent studies have considered more ternary oxides such as MgAl2O4 (Davies et al. 1994), La2CuO4 and Nd2CuO4 (Kenway et al. 1995) and where detailed experimental evidence is available, e.g. via Scanning Tunnelling Microscopy (STM) on WO3, the simulations reproduce the surface structure (Oliver et al. 1996). The scope of simulations has been extended further by considering surface microfaceting of MgO (Watson et al. 1996) and the energies of the physical and chemical adsorption of water onto the surfaces of MgO (de Leeuw et al. 1995, 1996). The calculation of surface structures and energies is important step toward modelling and understanding the behaviour of polycrystalline minerals and ultimately rock assemblages, because it is the interfaces that modify the behaviour of the polycrystalline material. The aim of the work described in this paper is to apply the simulation methods that have been successful in modelling simple ceramic oxide surfaces to a silicate mineral. This will allow us to investigate the extent to which potentials developed and used for simulating bulk silicate crystal properties can be used to model surfaces. The first silicate considered was forsterite (Mg2SiO4) because its bulk properties have been extensively modelled (Price et al. 1987a, b; Parker and Price 1989) and the THB1 potential is considered a reliable potential for predicting the bulk properties of this system. Surface methodology 1 The basic approach, used within the METADISE program (Watson et al. 1996), is to consider the crystal as a 71 The specific surface energy is defined as the energy per unit area required to transform a bulk region into a surface region. The surface energy is thus given by; g5 Usurf2Ubulk A where Usurf refers to the energy of region I of a surface calculation, Ubulk refers to the energy of an equivalent number of bulk ions, and A is the area. A condition for calculating the surface energy is that there is no dipole perpendicular to the surface in the repeat unit, because such a dipole leads to a divergent surface energy (Bertaut 1958). Tasker (1979a) defined three types of surfaces, I, II and III. Type I surfaces are composed of stoichiometric layers and thus have no dipole perpendicular to the surface (Fig. 1a). Type II surfaces are composed of multi-layer repeat units which have no dipole (Fig. 1b). Type III surfaces have a multilayer repeat unit and have a dipole perpendicular to the surface (Fig. 1c). Type III surfaces would therefore be predicted to be unstable, however recent work on NiO (Oliver et al. 1993), MgO (Watson et al. 1996) and MgAl2O4 (Davies et al. 1994) has shown that these surfaces can be stabilised by oxidation or reconstruction to remove the dipole. The calculated surface energies of a crystal are very difficult to compare with experiment, because of the lack of experimental data. However, one strategy for comparing the relative stability of the surfaces is to evaluate the morphology which can be either equilibrium or kinetic growth forms. Crystal morphology Fig. 1a–c Illustration of the three types of surfaces. a Type I, b Type II, c Type III Gibbs in 1878 (see Gibbs 1928) first proposed that the equilibrium form of a crystal should, for a given volume, possess minimal surface energy i.e. block consisting of a stack of planes periodic in two dimensions and parallel to the defect (surface). An assumption is made that only those ions close to the defect will need to relax from their lattice sites and thus the blocks are split into an inner region (region 1) which contains those ions close to the defect and an outer region (region 2) incorporating ions more distant from the defect. The ions in region 1 are allowed to relax to their mechanical equilibrium while those in region 2 are fixed at their bulk lattice sites. The energy and forces on the ions are calculated within the Born model of solids (Bom and Huang 1988) using a potential model which includes long range electrostatic terms and short range interactions. We have used the THB1 potential (Price et al. 1987a, b; Parker and Price 1989), which also includes oxygen polarisability via the shell model (Dick and Overhauser 1958) and directional bonding via three body forces (Sanders et al. 1984) and has been shown to reproduce successfully the structure, elastic properties and phonon frequencies of forsterite. g5Si gi Ai is a minimum for constant volume where gi and Ai are the surface energy and surface area of the ith crystallographic face. Wulff (1901), starting from the Gibbs proposal, suggested that the shape thus defined would be such that hi, the normal vector to the face from a point within the crystal, would be proportional to the surface energy of that face, gi; hi5l gi where l is a constant that depends on the absolute size of the crystal resulting in the crystal morphology as shown in Fig. 2. This expression will only hold for crystals grown with all faces in equilibrium and takes no account of kinetic factors such as growth rate. The difficulty is that it is often unclear when a crystal exhibits its equilibrium form. However, for small crystals where rearrangement of the crystal at each stage of the growth process is possi- 72 Fig. 2 Schematic representation of the calculation of crystal morphology ble, due to the small distances over which material has to travel, or if the crystal habit is formed over geological time scales, a morphology close to the equilibrium form might be expected. However, in the latter case the surface energy may be dominated by impurities, the concentration of which will be significant in geological materials. In an attempt to access kinetic factors we have also calculated the growth morphology. One phenomenological approach for estimating the growth morphology of crystals was introduced by Donnay and Harker (1937). In this, the morphological importance of a surface is inversely proportional to the interplanar repeat distance dhkl, i.e. hi5l 1 dhkl Thus low index surfaces with large interplanar spacing are predicted to dominate the morphology, which is in general true. A second approach is to use attachment energies, which are defined as the interaction of a layer of thickness dhkl with the rest of the crystal, assuming a bulk termination (Hartman and Bennema 1980). A large attachment energy indicates that ions coming onto the surface are strongly bound and so the attachment energy can be used as a measure of the growth rate of the surface. By finding the slowest growing cut for each surface, a pseudo-kinetic morphology can be predicted by relating the inverse of the surface’s morphological importance (i.e. the height of the normal) to the growth rate growth rate via the attachment energy. hi5l E(att)i Essentially the attachment energy approach is related to the Donnay-Harker scheme (Hartman and Bennema 1980). The use of the attachment energies adds atomistic information by adjusting the morphology to account for faster growth of those surfaces which have strong bonds perpendicular to the surface. A less favoured phenomenological approach is to argue that as the most stable surfaces are well known to grow the most slowly then the surface energies, which are a measure of the relative stability are proportional to the growth rates. In energetic terms the most stable surfaces prefer least new material attached to their surface. As the resulting morphology is equivalent to the equilibrium form this approach does provide a explanation for the good agreement between calculated equilibrium forms and large grown crystals. In practice these simple models can only be used as a guide because the growth rate will depend on many factors that are not included, for example, the free energy of formation of steps and kinks on each surface and the properties of the growth medium including saturation, pH and ionic strength. The low index surfaces of forsterite (Mg2SiO4) Forsterite (Mg2SiO4) is an end member of the mineral olivine, which shows a complete solid solution from forsterite, the magnesium end member, to fayalite (Fe2SiO4), the iron end member, with the typical sample of the mineral olivine (Mg, Fe)2SiO4 containing approximately 10% Fe (Fo90). It has orthorhombic symmetry (space group Pnma with a510.20 Å, b55.98 Å, c54.75 Å), and consists of independent SiO4 tetrahedra linked by divalent cations in octahedral co-ordination. The divalent cations occupy two distinct crystallographic sites with the oxygen sub-lattice deviating significantly from the ideal close packing arrangement, resulting in variation in Mg-O bond distance (2.05–2.21 Å). The large deviation of the oxygen sub-lattice results in forsterite having a very complex set of crystal planes. For the simplest surfaces, those containing plane indices of either 0 or 1 (termed the 1 index surfaces), there are seven distinct sets of planes ({100}, {010}, {001}, {110}, {101}, {011}, {111}) because of the crystal’s orthorhombic symmetry. Figure 3 shows the first non dipolar repeats found for each surface with indications of addition planes which may be cut to form non dipolar repeats and include both type II and type III (or dipolar surfaces). The type III surfaces can in all cases be reconstructed with a zero dipole by moving half of the ions from the top to the bottom of the repeat (e.g. {010a}). The repeat units were generated assuming that the tetrahedra would remain intact, as previous calculations (Davies 1992) have indicated that the breaking of the Si-O interaction always results in less stable surfaces. This causes some of the cuts not to be simple cleavage planes as the oxygens are cut according to the position of the silicon they are bonded to. For example, the {111a} surfaces show an O1 and O2 above the Mg2 which was on the plane of the cut due to them forming part of the tetrahedron of the silicon below this magnesium. Some of the surface cuts are complicated by cation arrangement at the surface. For example, the {101b} surfaces cut a plane of four magnesium ions in half resulting in 50% occupancy of the surface magnesium layer. Therefore there are six possible cation arrangements at the surface, but in general this can be reduced, in this case to three, by considering surface symmetry. 73 Fig. 3 Schematic representation of the surface stacking sequences for the 1 index surfaces of forsterite, (100, 010, 001, 110, 101, 011, 111). The figure shows the sequence of ions perpendicular to the the surface, and the possible cuts for each surface that result in a zero dipole repeat unit A general problem of evaluating the surface energy is ensuring its convergence and is achieved by simulating a sufficiently large crystal region. Table 1 shows the variation in the relaxation displacement of the first 28 ions as a function of depth for the {100a} and {010b} surfaces with region 1s containing 112 ions (176 species). This table illustrates two points. The first is that the relaxation is surface-dependent and continues over 6 Å into region 1. The second is that because of forsterite’s orthorhombic nature, the area of the surface unit cells vary according to the orientation and therefore the same number of ions results in varying region 1 depths for different surfaces. This combination means that great care must be taken to ensure that enough ions are considered to ensure 74 Table 1 The displacement on relaxation (x, y, z in Å) of the top 28 cores from their bulk positions in the METADISE minimised (region 15112 ions) 100a and 010b surfaces as a function of depth Table 2 Lowest unrelaxed and relaxed surface energies for the 1 index surfaces of forsterite a The 110b surface has a relaxed surface energy of 1.92 Jm22. b The 3 indicates this to be the third of three possible cation arrangements. c The 111a surface has a relaxed surface energy of 1.90 Jm22. 100a Ion Mg2 O1 O1 Si O2 O3 Mg1 Mg1 O3 O2 Si O1 O1 Mg2 Mg2 O1 O1 Si O2 O3 Mg1 Mg1 O3 O2 Si O1 O1 Mg2 010b Depth Displacement (Å) (Å) X 0.00 0.54 0.54 1.22 1.28 1.74 2.22 2.22 2.71 3.17 3.23 3.91 3.91 4.45 5.12 5.66 5.66 6.34 6.40 6.86 7.35 7.35 7.83 8.29 8.35 9.03 9.03 9.57 0.28 20.04 20.04 20.06 0.15 0.05 0.02 0.02 20.05 20.04 0.05 0.05 0.05 20.04 0.02 20.02 20.02 20.01 0.00 0.00 20.01 20.01 20.01 20.01 0.00 0.00 0.00 0.00 Y Ion Z 0.00 0.06 20.06 0.00 0.00 0.00 20.04 20.04 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 20.09 0.11 0.11 20.02 20.06 0.07 20.03 20.03 20.07 0.02 0.00 20.01 20.01 20.02 0.00 20.01 20.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 Mg1 O1 O1 O3 Mg2 O2 O3 Si Mg2 Si O2 O1 O1 Mg1 Mg1 O1 O1 O2 Si Mg2 O2 O3 Si Mg2 O3 O1 O1 Mg1 Depth Displacement (Å) (Å) X 0.00 0.21 0.21 1.50 1.50 1.50 1.50 1.50 1.50 1.50 1.50 2.79 2.79 2.99 2.99 3.20 3.20 4.49 4.49 4.49 4.49 4.49 4.49 4.49 4.49 5.78 5.78 5.99 0.39 0.36 0.51 20.63 0.29 0.14 20.17 0.01 20.16 0.00 20.09 0.07 0.24 0.28 20.24 20.01 20.07 0.07 0.06 0.00 20.15 20.23 20.16 0.13 0.31 20.03 20.10 20.15 Y 0.58 20.98 20.69 20.17 20.15 20.44 20.38 20.39 20.26 20.12 20.03 20.14 0.50 20.51 20.17 20.36 0.22 20.15 20.22 20.11 20.05 0.00 20.04 20.09 20.16 20.02 20.11 20.11 Z 20.02 0.60 20.40 20.47 0.23 0.25 0.43 0.29 20.15 20.23 20.28 0.11 20.37 0.04 0.20 0.05 20.17 0.04 0.04 0.09 0.01 0.08 0.03 20.01 20.04 0.05 0.01 0.02 Plane Index Miller Index Unrelaxed surface energy (Jm22) (cut) Relaxed surface energy (Jm22) (cut) Attachment energy (eVySi) 100 010 001 110 101 011 111 200 020 002 220 101 011 111 2.23 3.84 4.61 6.56 3.51 4.88 4.01 1.28 1.61 2.02 1.90 1.70 1.81 1.80 21.98 27.20 28.58 214.06 27.02 25.47 25.21 convergence of the surface energy. This is further illustrated by the surface energy of the {100a} surface which converges for a region size of 1 unit cell (44 species, 10.24 Å) while the {101b} requires 3 unit cells (132 species, 13 Å) and the {010a} requires 4 unit cells (176 species, 23.95 Å). Minimisation of the 1 index surfaces gave rise to the lowest surface energies for each index before and after minimisation shown in Table 2. The relaxation results in only a small change in the order of stability of the surfaces with the order before minimisation of, 100 ,101 ,010 ,111,001,011,110 and after minimisation, 100 ,010,101,111©011,110,001. (a) (b) (b) (a) (b3)b (a) (b) (a) (b) (a) (a)a (b3)b (a) (b)c The 2 index surfaces On inclusion of surfaces with plane indices of up to 2, the total number of structurally unique orientations increases from 7 to 19. The additional 12 have much larger surface areas, and thus for a given number of ions produce surface unit cells of less depth. This greatly increases the computational effort for three reasons. (i) The number of ions required to ensure convergence of the surface energy increased. For example, the {012} surfaces for a cell of 44 species results in a depth perpendicular to the surface of only 2.54 Å requiring 308 species (a depth of 17.8 Å) to ensure convergence in METADISE. (ii) The occurrence of many type III cuts which on reconstruction to form a zero dipole repeat results in 50% 75 occupancy of the surface cation sites and hence a number of possible cation rearrangements need to be considered. (iii) The occurrence of asymmetric stacking sequences. This occurs when a zero dipole repeat unit is not symmetric which results in different surfaces at each end of the repeat unit. Fortunately for forsterite where such stacking sequences are found they are composed of two non-dipolar units alternating and thus a mirror image of the stacking sequence can be formed (termed asymmetric inverse). Thus allows growth of both surfaces cut at both ends of the crystal, and thus the surfaces can be treated independently (Fig. 4). Minimisation of all the 2 index surfaces was achieved although some of the more unstable surfaces showed a number of high energy minima which was taken to be an indication of the surface instability. For example, the {210a} before minimisation have a surface energy of 7.5 Jm22 and relaxed to surface energies at both 2.39 Jm22 and 2.60 Jm22. This type of multiple minimum was investigated by varying the number of steps between updating the second derivative matrix and adjusting the maximum displacement an ion can move in a single step. In general the more stable surfaces achieve minimisation to a global minimum while some of the unstable surfaces showed multiple minima. In each case we believe the lowest energy surface has been identified. The lowest surface energies for each index are shown in Table 3 indicating that in general the symmetric surfaces are more stable than the asymmetric. Crystal morphology The order of surface stability including all of the one and two index surfaces changes from 100,201,101,121,021,221,010,111,211,210 ,212,011,122,110,001,012,120,112,102 before relaxation to 100,201,221,211,010,101,111,011,121,021 ,120,110,210,001,212,122,122,012,102 Fig. 4a–d Illustration of asymmetric and asymmetric inverse stacking sequences for the {221}b surface. a 221bas repeat unit, b 221basi repeat unit c and d the two alternating non dipolar repeats within the surface repeat Table 3 Lowest unrelaxed and relaxed surface energies for the 2 index surfaces of forsterite after relaxation. This change in stability has an effect on the calculated equilibrium morphology. Figure 5 shows a comparison of the Donnay Harker (Fig. 5a), attachment energy (Fig. 5b) and equilibrium morphologies, before and after relaxation (Fig. 5c, 5d), with two observed morphologies for natural olivine (Figs. 5e, 5f) (Dana and Ford 1958; Deer et al. 1982). Plane Indexy Miller Index Unrelaxed surface energy (Jm22) (cut) Relaxed surface energy (Jm22) (cut) Attachment energy (eVySi) 012 021 102 112 120 121 122 201 210 211 212 221 7.00 3.70 8.15 5.85 7.18 3.69 5.11 2.56 4.67 4.29 4.75 3.72 2.29 1.83 2.30 1.99 1.89 1.82 2.12 1.56 1.95 1.59 2.02 1.58 216.20 27.41 231.93 216.06 215.78 27.75 216.03 23.38 25.69 26.87 29.88 28.30 (b) (a) (c) (a) (b) (a) (c) (a) (b) (a) (b) (b1) (a) (a) (a) (a) (b) (a) (c) (a) (b) (a) (b) (b1) 76 absent. In addition the {210} surfaces are also absent which has a large impact on the observed morphologies. The unrelaxed equilibrium morphology (Fig. 5c) shows a number of features consistent with those observed. The {100}, {201}, {021} and {010} are all present, however the {121} and {221}, which are not observed, also appear. This time the Donnay Harker and attachment energy schemes the {201} is present however the {101} is now absent. The relaxed equilibrium morphology expresses the {201} in agreement with observations unlike the Donnay Harker and attachment energy and the {101} unlike the unrelaxed equilibrium morphologies and indicates the importance of surface relaxation on surface stability. It does however, include a number of additional surfaces, namely the {221} and {211} in contradiction to the observed morphologies. In summary, none of the morphologies match the observed samples in all respects. The Donnay-Harker scheme is the only one to match those observed with respect to the appearance of the {210} surface while the relaxed equilibrium morphology is the only scheme with correctly predicts the appearance of both the {201} and {101} indicating the relaxation is an important factor. Surface relaxation and structure Fig. 5a–f Comparison of calculated a Donnay Harker, b attachment energy, c unrelaxed, d relaxed and observed, e Dana (1958), and f modified after Deer et al. (1982) crystal morphologies of forsterite Matching the predicted morphology with those grown in the laboratory can be difficult despite the ideal conditions. Kinetic factors play a major role in influencing the morphology because of fast growth rates. However, even for such crystals thermodynamics may play a part especially if the crystals selected for study are those that keep the same morphology with time. For mineralogical samples the factors controlling crystal growth are much harder to determine. Additionally there are increase amounts of impurities which can segregate and hence modify the growth characteristics of different surfaces to different extents. Comparison of the different models shows that the Donnay-Harker morphology (Fig. 5a) gives good agreement with those observed (Figs. 5e, 5f) and has similarities with the relaxed equilibrium morphology (Fig. 5d). It expresses all of the required surfaces except the {201}, which is absent because of the dominance of the low index {101} surface causing the crystal to appear tabular. The attachment energy morphology (Fig. 5b) shows some agreement with observed morphologies although the agreement is not as good as that for the Donnay Harker scheme. The {100} surfaces are too dominant leading to a plate like crystal, with the {201} surface once again Table 2 clearly shows that the {100a} surfaces are the most stable before (2.23 Jm22) and after (1.28 Jm22) relaxation. The small relaxation energy implies that the surface cut does not remove a large number of bonds and this is seen in Fig. 6a, which shows the {100a} surface before and after relaxation. The high surface density with alternate SiO2 and MgO layers terminating in MgO gives rise to a surface structure and energy very similar to that of the {100} surface of MgO (1.20 Jm22). Relaxation is very minor with the surface magnesium layer relaxing slightly (0.3 Å) towards the bulk, with the second layer (SiO2) relaxing by a very small amount (,0.1 Å) away from the bulk, resulting in a slight undulating surface topology. This is in agreement with the only experimental study of the surface structure of olivine (Hochella 1990), which looked at the {100} surface using Low Energy Electron Diffraction (LEED). They observed a 131 surface unit cell with only small relaxations. The next three lowest relaxed surface energies, the {201a}, {221bas1} and {211a} are all within 0.03 Jm22 with the {201} surfaces the most stable before and after minimisation and is the result of, like the {100a}, a compact surface with magnesium in the top layer and the first silicon 1.2 Å below the surface. This results in a very stable surface with only minor relaxations of the surface magnesium and oxygens (Fig. 6b) and gives rise to a definite stepped surface topology. The {221bas1} is an asymmetric surface which initially has a single oxygen pointing straight up as a result of keeping the tetrahedra complete. The initial surface en- 77 Fig. 6a–f The structure on the three lowest energy surfaces of forsterite before and after relaxation. a 100a, b 201b, c 221bnsl, d 211b, e 010b and f 101b3 ergy is surprisingly low (3.72 Jm22) and is a result of the rest of the surface displaying a generally flat termination with predominantly magnesium rather than silicon at the surface. Relaxation of this surface is quite dramatic with surface tetrahedra rotating to allow the isolated oxygen ion to co-ordinate to two surface magnesium atoms (Fig. 6c). The surface relaxes to accommodate this with the surface oxygens moving as much as 2.00 Å resulting in an undulating surface structure. The next surface, the {211a}, has a much less compact jagged structure with oxygen ions sticking directly out from the bulk. On relaxation one magnesium ion moves across the surface by 1.5 Å to coordinate to an isolated oxygen (Fig. 6d). This appears to distort the structure but only results in slight rotation of some of the surface SiO4 tetrahedra with the oxygen ions moving to increase their interactions. Relaxation of up to 0.1 Å continues for up to 7 Å into the surface and the whole result is a large stabilisation of the surface from 4.29 Jm22 to 1.59 Jm22. Figure 6e shows the surface structure of the {010b} which is a type III surface formed by cutting a layer of two magnesium ions in half. This results in a surface with an isolated magnesium ion which on relaxation moves both toward the bulk by 0.4 Å and laterally within the plane of the surface by 0.6 Å moving away from one oxygen and toward another. This type of lateral move- 78 ment is observed up to 6 Å into the surface and together with rotation of the SiO4 tetrahedra gives rise to a flatter and more compact, but distorted, surface structure. The {101b3} surface (Fig. 6f) is the result of cleaving a plane of four magnesium ions. This leads to three independent surfaces for which the surface energy varies from 7.00 Jm22 to 3.51 Jm22 for the unrelaxed and 2.28 Jm22 to 1.70 Jm22 for the relaxed surfaces. The figure shows that the initial cut is quite jagged giving rise to exposed magnesium ions at the surface. On relaxation these magnesium ions move towards the bulk occupying positions between oxygen ions increasing the number of nearest neighbour oxygens. This again results in an almost flat surface structure which is consequently relatively stable. Conclusions This aim of this paper was to demonstrate the scope and limits in using computer simulation to model the surface energies and structures of a silicate mineral. The results of these simulations have provided predictions of the surface structures which we hope will be verified by experimental study. Indeed, the recent success of atomic force microscopy (AFM) on a range of mineral surfaces suggests this is now possible and that with reliable models of the surface from simulation studies interpretation of these images may be easier. The basis for suggesting we have a reliable model is first, that the potential model chosen has previously been shown to predict accurately a wide range of bulk crystal properties (Price and Parker 1987a, b; Parker and Price 1989). Secondly, there is good accord with the crystal morphology of natural olivine, although for reasons discussed above an exact match is not yet possible because of the limitations of the kinetic growth models. Our future work in this area will include the application and development of more sophisticated growth models allowing not just qualitative comparison with morphology but predictions of the influence of surface defects and impurities. The latter will be particularly important when understanding the defect properties of geological materials. Acknowledgements We would like to thank the EPSRC and the NERC for financial support and Molecular Simulations Inc. for the provision of INSIGHT II. References Bertaut F (1958) The electrostatic term of the surface energy. CR Acad Sci 246: 3447 Bom M, Huang K (1988) Dynamical theory of crystal lattices. Clarendon Press, Oxford Colbourn EA, Mackrodt WC, Tasker PW (1983) The segragation of calcium ions at the surface of magnesium oxide: theory and calculation. J Mater Sci 18: 1917–1924 Dana ES, Ford WE (1958) A textbook in mineralogy. Wiley, New York Davies MJ (1992) Computer simulation of oxide surfaces and grain boundaries. PhD thesis, University of Bath Davies MJ, Parker SC, Watson GW (1994) Atomistic simulation of the surface structure of spinel. J Mater Chem 4: 813–816 Deer WA, Howie RA, Zussman J (1992) An introduction to the rock forming minerals. Longman, London de Leeuw NH, Watson GW, Parker SC (1995) Atomistic simulation of the effect of dissociative adsorption of water on the surface structure and stability of calcium and magnesium oxide. J Phys Chem 99: 17219–17225 de Leeuw NH, Watson GW, Parker SC (1996) Atomistic simulation of adsorption of water on three-, four-, and five-coordinated surface site of magnesium oxide. J Chem Soc Faraday Trans 92: 2081–2091 Dick BG, Overhauser AW (1958) Theory of the dielectric constants of alkali halide crystals. Phys Rev 112: 90–103 Donnay JD, Harker G (1937) A new law of crystal morphology extending the law of Bravais. Am Mineral 22: 446–467 Gibbs JW (1928) Collected works, Longman, New York Hartman P, Bennema P (1980) The attachment energy as a habit controlling factor. J Cryst Growth 49: 145–156 Hochella MF Jr (1990) Atomic structure, microtopography, composition and reactivity of mineral surfaces. In: Ribbe PH (ed) Mineral-water interface geochemistry (Reviews in Mineralogy vol 23). Mineralogical Society of America, Washington DC, pp 87–132 Kenway PR, Parker SC, Mackrodt WC (1995) The calculated structure, stability and composition of the low index surfaces La2CuO4 and Nd2CuO4. Surf Sci 326: 301–310 Mackrodt WC, Stewart RF (1977) Defect properties in ionic solids: I point defects at the surfaces of face centred cubic crystals. J Phys C 10: 1431–1446 Oliver PM, Parker SC, Mackrodt WC (1993) Computer simulation of the crystal morphology of NiO. Modelling in Materials Science and Engineering 1: 755–760 Oliver PM, Parker SC, Edgell RG, Jones FH (1996) Computer simulation of the surface structures of WO3. J Chem Soc Faraday Trans 92: 2049–2056 Parker SC, Price GD (1989) Computer modelling of phase transitions in minerals. Advances in Solid State Chemistry 1: 295– 327 Price GD, Parker SC, Leslie M (1987a) The lattice dynamics of forsterite. Mineral Mag 51: 157–170 Price GD, Parker SC, Leslie M (1987b) The lattice dynamics and thermodynamics of Mg2SiO4 polymorphs. Phys Chem Minerals 15: 181–190 Sanders MJ, Leslie M, Catlow CRA (1984) Interatomic potentials for SiO2. J Chem Soc Chem Commun:1271–1273 Tasker PW (1978) U.K.A.E.A Report, AERE-R-9130, United Kingdom Atomic Energy Authority, Harwell, UK Tasker PW (1979a) The stability of ionic crystals. J Phys C 12: 4977–4983 Tasker PW (1979b) The surface energies, surface tensions and surface structure of the alkali halide crystals. Philo Mag 39: 119–136 Tasker PW, Duffy DM (1984) The structure and properties of the stepped surfaces of MgO and NiO. Surf Sci 137: 91–102 Watson GW, Kelsey ET, de Leeuw, Harris DJ, Parker SC (1996) Atomistic simulation of dislocations, surfaces and interfaces in MgO. J Chem Soc Faraday Trans 92: 433–438 Wulff G (1901) Zur Frage der Geschwindigkeit des Wachstums und der Auflosung der Krystallflachen. Z Kristallogr Kristallgeom 34: 949