Furuday Discuss., 1993,95, 75-84 Computer Modelling of Inorganic Solids and Surfaces Stephen C. Parker,* E. Toby Kelsey, Peter M. Oliver and James 0.Titiloye University of Bath, School of Chemistry, Claverton Down,Bath, UK BA2 7AY We have recently performed calculations on a number of inorganic solids and their surfaces. The emphasis of this work is to understand the role of surface defects, either intrinsic defects or additives, in modifying structure and stability at the atomic scale. The basis of the approach is to use energy minimisation to obtain the most stable configuration. In this paper we describe simulations on NiO which predict that the { 11 l} configuration is stabilised by surface oxidation at high temperatures. Further simulations of calcite model the effect of additives on morphology, and we describe the effects of lithium, magnesium and hydrogenphosphate _additives.We find that both magnesium and phosphate stabilise the { lOlO} surface while lithium stabilises the basal plane. Finally, we present preliminary work on calcite and barium sulfate which implies that these methods may provide useful insights in nucleation and crystal growth at high ionic strengths. A detailed description of surface structure and reactivity is central to understanding crystal growth at the atomic scale. This is particularly important for interpreting the role of additives which is of concern for many technological processes. However, experimental studies of the structure and stability of interfaces, particularly in polar solids, are often difficult.Atomistic simulation provides a complementary technique for probing the surface structure, stability and behaviour of defects. Atomistic simulation is now becoming well established and stems from the develop- ments of Tasker and Mackrodt, who demonstrated the potential of this technique by successfully modelling the surface properties of a range of ceramic oxides which include the cubic rock-salt oxides MgO, CaO and NiO19* and more recently materials such as L~,CUO,.~However, it is only very recently that these techniques have attempted to rationalise the effect of additives to the growth and morphology of inorganic solid^.^ The major complication with most inorganic solids is their ionic character which has a profound influence on the stability of their surfaces and hence on their morphology.However, the development of the techniques and increase in available computer power now make simulation of polar solids possible. We will attempt to highlight this by describing our current work on NiO and CaCO,. Experiments on each material show that defects and impurities significantly modify the observed morphology. Hence we aim to discover whether these simulation techniques can model this observed behaviour. In particular, we need to consider whether simulations, based on calculations of free surfaces, can be applied to the morphology of a growing crystal at either high temperatures or in the presence of a solvent containing varying amounts of impurities. However, before discussing the results we will briefly summarise the simulation approach.Theoretical Methodology The atomistic simulations are based on the Born model of ionic solids, in which the ions interact via long-range electrostatic forces and short-range forces. The latter include a representation of the repulsion between neighbouring charge clouds, and a shell-model 75 Computer Modelling of Surfaces Table 1 Interatomic potential parameters: Buckingham potential: V = A,, exp( -r,,/p,,)-CJr6 (i) calcite (shell-ion model) Ca-0 1111.0 0.30926 0.0 Ca: 2 + Mg-0 583.754 0.29654 0.0 Mg: 2 + Li-0 5 10.092 0.29654 0.0 Li: 1 +o-*.o 6959.0 0.23654 0.0 0:-0.97310c-0 527.5 0.1550 0.0 C: 0.91930 three-body potential V = 112 K(0-Oo)2 0-C-0 K= 2.6, 0, = 120.0O oxygen core-shell spring constant 74.92 (ii) barium sulfate (rigid-ion model) Ba-0 7977.5 0.2698 0.0 Ba: 2 +o.-*o 360 10.0 0.19756 0.0 0: -1.0s-0 1827.9 0.19910 0.0 s:2+0 0-S-0 K = 9.09724, 0, = 109.47O (I Ref.24. description of the ions to account for their p~larisation.~ However, for complex materials, particularly those containing polyanions, it is often necessary to include angle-dependent forces to account for the covalence, as is the case for calcium carbonate.6 All of these forces, collectively called the potential model, are specified by simple parametrised analytical functions. The reliability of the simulations are dependent on the choice of parameters. However, sometimes the reliability of the parameters for a particular compound is compromised to ensure transferability between a number of systems.In all cases considerable effort must be made to ensure that the parameters are sufficiently reliable for the particular application. The parameters for nickel oxide and calcite were obtained by fitting directly to experimental data, such as the structural, elastic and dielectric data. The parameters for the preliminary study using barium sulfate were fitted to a range of sulfates, for the latest parameters see Allan and Rohl (this Discussion). The parameters we used for NiO are those of Sangster and St~neham,~ and for calcite and barium sulfate are given in Table 1. The low energy and hence the most common surfaces of a crystal are generally those of low Miller index.These planes are the closest packed with large interplanar spacings. However, in ionic crystals other constraints also apply. If the Madelung sums are not to diverge with increasing crystal size, then the crystal must not only be electrically neutral but also have no net dipole moment perpendicular to the surface. Bertaut’ demonstrated that when there is a dipole moment perpendicular to the surface, the surface energy diverges and is infinite. Such surfaces are therefore unstable. This result provides a simple method for selecting the potentially important crystal faces for simulation. Thus the simulation strategy is first to cleave a crystal to generate a slab or block which is charge neutral and non-dipolar perpendicular to the surface. The crystal faces fall into three categories, see Fig.1. Type I surfaces comprise planes containing all component ions in their stoichiometric ratio. Common examples include the (100) and (1 lo} surfaces of rock salt which have equal numbers of cations and anions. Thus cleaving the crystal at any position will yield a non-dipolar repeat. Type I1 surfaces require more care as these must be cut at a specific plane to avoid a dipole. In the example in Fig. 1(b),the cut must produce an X-M-X repeat rather than M-X-X or X-X-M. Type 111 surfaces if cut at any plane produce a dipole moment perpendicular to the surface, and cannot occur naturally without the adsorption of foreign atoms or surface roughening.As a consequence they are usually S. C. Parker et al. .o.o.o.o] neutral non-dipolar repeat unit (b).... ~o~o~o.o 0 M2+....0.-(c) dipolar repeat unit Fig. 1 Schematic representation of (a) type I, (b) type I1 and (c) type I11 surfaces the least stable.' One apparent exception is NiO in which the morphology is dominated by the (1 11) surface when grown at elevated temperatures. The final structure (and energy) of a surface is then determined by the requirement that the system is in mechanical equilibrium. This is achieved by allowing the ions to relax to the point at which they experience zero net force. The number of relaxed surface planes is increased until the surface energy converges. Defects and impurities can be accommodated in the surface but only such that the net charge is zero.The equilibrium crystal morphology of a material can be determined by applying Wulff's Theorem,8 although it was Gibbs9 who first proposed that the equilibrium form of a crystal should possess minimal total surface energy for a given volume. Wulff proposed that a polar plot of surface energy as a function of the orientation of normal vectors would give the crystal morphology. This assumes that the crystal is small and can rearrange during growth due to the short distances which ions have to travel. Thus on calculating the surface energies we can determine an equilibrium morphology which can be compared with experiment. Results and Discussion Surface Energies and Morphology of NiO We are attempting to address the ambiguity whereby NiO expresses the notionally unstable (1 1 1) surface (type 111) in collaboration with W.C. Mackrodt'O by examining the surface energies of NiO and comparing the resulting morphology with that determined by Handwerker et aZ.l experimentally. The energies of the three lowest-index surfaces were calculated. As noted above, the {loo}and { 1 lo} surfaces can be created unambiguously by cleaving the rock-salt structured NiO as they are type I surfaces. In contrast the simple cleavage of the { 11 l} results in a dipole perpendicular to the surface, i.e. a type I11 surface. The reconstruction for stabilising the ( 1 1l}plane which maintained stoichiometry involved removing one half of the surface atoms, for example anions and adding them to the opposite cation surface. On calculating the surface energies for the { 1001, { 1101, and { 1 11) planes the results agreed with the earlier work of Tasker2 and showed that as with MgO the most stable surface by far is {loo}.Hence the resulting morphology would be predicted to be cubic. However, one major difference between NiO and MgO is that nickel oxide is known to be Computer Modelling of Surfaces non-stoichiometric. Thus we can consider a reaction in which the surface reacts with oxygen to form holes (with a charge of 1 +) and charge-compensating defects such as nickel vacancies, with an effective charge of 2 -. Thus the reaction becomes: Niki + to2-+ V& + 2h’ + NiO (1) in Kroger-Vink notation Ni& represents a Ni2+ lattice ion, V!i a doubly charged Ni vacancy and h‘ an electronic hole.There is still considerable controversy relating to the character of the hole. The most stable electronic configuration from the atomistic simulations was a hole localised on the nickel, i.e. Ni3+. The total energy of the above reaction in the bulk was calculated to be 1.6 eV, which agrees well with previous work.I2 Thus we applied the reaction to the NiO surfaces which formed neutral clusters consisting of two holes and one vacancy. The surface energies as a function of the vacancy-hole cluster are given in Fig. 2. The most obvious feature is that the variation of surface energy with defect concentration is significantly different for each surface. The { 100) plane becomes less stable, the ( 110) plane has a minimum at 50% coverage, while the { 1 1l} plane shows the most complex behaviour with a local minimum at 50% and the lowest energy at 100% coverage.However, this distinctly non-Langmuir behaviour is not particularly unusual and has been predicted for a number of systems, most notably segregation of alkaline-earth-metal cations to the surfaces of magnesia.’ The results can be explained, in part, by steric factors as (100) is the most densely packed surface and hence the addition of defects serves only to destabilise the surface, while ( 1 10) is less dense and hence can accommodate a defect concentration of 50% before steric repulsion destabilises the surface. The (1 1l} plane has the lowest density of surface ions and hence the surface energy is lowered to the greatest extent.The local minimum at 50% corresponds to the complete removal of the half-plane of surface cations which were left after the necessary reconstruction required to remove the dipole. We next calculated the morphology using the surface energies for a given surface oxidation, and compared the resulting morphologies of Handwerker et al.” We found good agreement at 75% coverage, Fig. 3. Thus we interpret the observed morphologies as being bulk NiO with a surface hole concentration corresponding to the 75% defect coverage. It is important to emphasise that a high surface-defect concentration does not imply a high bulk concentration, as even a few ppm concentration in the bulk can segregate and fully cover the surfaces.c 0.5-I I I I I I I I I0.0 0.0 12.5 25.0 37.5 50.0 62.5 75.0 87.5 100.0 coverage Fig. 2 Surface energy vs. coverage for the (100) (O),{l10) (D) and ( 111) (A)surfaces of NiO S. C. Parker et al. Fig. 3 (a) Experimental morphology, (b) calculated morphology at 75% coverage for NiO Impurity Segregation to the Surfaces of Calcite Much of our current effort is directed towards modelling the effect of additives on surface structure and morphology. One chosen example was calcite because there is considerable experimental evidence with which to test the reliability of our approach13-15 and because there has been considerable success of Mann's group at Bath in revealing insights into the nature of calcite surfaces.I4 In this section we describe simulations modelling the influence of Mg2+, Li+ and HPOi-additives on the surface stability and morphology.The calcite structure is related to the rock-salt structure and can be envisaged as a replacement of the Na atoms by Ca and the Cl atoms by CO, groups coupled with a compression of the rhombohedron along the three-fold axis. Experimental work on calcite has shown that the low-index surfaces dominate the crystal morphology.16 A schematic representation of the low-index surfaces of calcite is shown in Fig. 4. The surface energies for the low-index faces { 1Oi4},{ lOiO}, { 1120} and (0001) were calculated using the energy- minimization method described above are given in Table 2.These results show that the { 1Oi4} surface, which has the smallest surface area, is the most stable face, while the first order prismatic faces { lOiO} and { 1 lzO} are about equal in stability. Of the low-index faces the (0001) face was found to be the least stable. This is not surprising as the basal plane is comprised of alternate layers of Ca2+ and C0:-ions in successive planes, revealing it to be a type I11 surface. The most stable reconstruction was similar to that of the NiO { 11 l} in which the surfaces were terminated with a half plane of cations Fig. 4(d). Another interesting feature of this system is that there was little surface relaxation particularly for the { 1074) surface. This indicates that the pure surfaces were simple bulk terminations.This was confirmed by the rhombohedra1 morphology of the { 10741 surface predicted from the surface energies which agrees well the experimentally observed morphology, Fig. 5(a).17We next added impurities to the surfaces, first magnesium and lithium cations and secondly, the HPOi-anion. (c) {lilO] t Fig. 4 Schematic representation of the stacking sequences for the low-index planes of calcite Computer Modelling of Surfaces Table 2 Calculated energies for selected low-index surfaces of calcite with and without defect substitution defect surface segregation energy energy /kJ mol-surface surface energy (pure) /J m-* /J m-? Mg2+ Li + Mg2+ Li + { 1oi4) { 1120) (ioio} (0001) 0.23 0.50 0.52 1.60 1.28 1.37 0.84 1.14 1.28 1.08 0.76 0.32 127.3 130.2 84.9 115.7 128.3 87.7 65.6 -208.3 The addition of Mg2+ or Li+ ions at the surfaces was carried out by substitution of the Ca2+ ions at the surface and comparing with substitution in the bulk.This is the segregation energy. The heats of segregation of Mg2+ at 100% coverage of different surfaces are shown in Table 2. A positive segregation energy indicates a preferential dissolution of magnesium into the bulk calcite crystal lattice. The results suggest that magnesium ions prefer to dissolve into the bulk from the stable faces. This result is not surprising given that magnesium and calcium carbonates can form a solid solution. However, unlike the high- temperature conditions in which NiO was grown in the previous example calcite is grown in a comparatively low-temperature solution.Thus given a high concentration of magnesium at the surface there will be little diffusion into the solid, and hence magnesium will remain at the surface. Furthermore, we have some confidence that a high concentration of magnesium ions will be present at the surface because calculations of the stability of magnesium and calcium in solution compared to the solid surface suggest that magnesium will segregate from solution to the surface, and unlike in the solid there should not be a large barrier to diffusion to the solid surface. The simulations predicted that incorporation of Fig. 5 Predicted crystal morphologies showing (a) { 10T4) rhombohedral; (b) { 1010) face stabilised with Mg2+;(c) (0001) face stabilised with Li+ and (6){ lOT0) face stabilised with HPOa- additives S.C. Parker et al. Mg2+ stabilises the first-order prismatic faces (lOTO} relative to other faces (Table 2). Thus the crystal morphology was modified to give a first-order prism capped with rhombohedra1 end faces as shown in Fig. 5(b).This morphology has similar features to that observed experimentally. 16~18 The substitution of Ca2+ ions with Li+ ions results in an effective negative charge which was compensated by the addition of a Li+ interstitial. We also considered substituting calcium with Li+ at two interstitial sites, but found in each case the lithium interstitial was destabilising. This is further exemplified by the large segregation energies from the calcite bulk to all surfaces, Table 2, by showing that replacing calcium by lithium is highly unfavourable. The calculated defect surface energies (Table 2) showed that all surfaces except the polar (0001) surface were destabilised.On cleaving the (0001) surface of pure calcite there are two coplanar Ca2+ sites where only one site is occupied. Hence by removing Ca2+, we can occupy both cation sites with Li+. This has the added effect of increasing the cation density at the surface which effectively stabilises the surface. The morphological consequence of the increased stability of the (0001) face is to predict a tabular crystal habit comprising basal (0001) and (1014) side faces, see Fig. 5(c), which agrees with recent experimental studies of calcite crystallization in the presence of Li+ ions.19 We next used the same procedure for simulating the effect of the hydrogenphosphate anion, namely substituting the surface carbonates by HPOi-, the potential model is given in Table 3.The only difference was that 100% replacement was not energetically feasible due to the size of the HPOZ- anion, but was with 50% coverage of each surface. The HPOi- anion was observed to prefer the { 10lo} surface energetically compared to ( 1 120) or { 1014). Fig. 4(4 shows the morphology predicted for the effect of HPOf- anions on calcite surfaces. The morphology is that of a first-order prism face, { lolo}, similar to that predicted for the Mg2+ additives. This leads us to question whether a monovalent anion will have the same effect as lithium and stabilise the (0001) surface.A further unexpected result on simulating the addition of hydrogenphosphate anions was the effect on the structure of thin films which, if it is a general effect, may enable simulation to be used to provide a coarse screening of potential additives used to inhibit nucleation. Applications to Nucleation Studies Simulations on thin films of ca. 25 A thick expressing the { 1070)surface (Fig. 6) show that adding hydrogenphosphate ions to the surfaces causes major disru tion to the crystal surface. However, on increasing the thickness of the film to ca. 75 K we find that after relaxation there is no disruption of the lattice. These results suggest that the hydrogenphos- phate ion inhibits CaCO, nucleation by forcing apart the molecules in the bulk and Table 3 Interatomic potential used for the hydrogenphosphate additives on calcite (i) Buckingham potential: V = A, exp( -r,/p,) -C,,/r6o.*.o 1388.773 0.362 3 1 175.0 25 P-0 9034.208 0.19264 19.8793 25 0-0, 3627.500 0.28749 3.47 this work 0-Ca 2396.800 0.29202 0.0 this work 0-H (ii) Morse potential: V &, = 5.896 eV, a = 26 Charges on P = + 3.4,O = -1.46667 and OH = -1.O.0,is the carbonate oxygen. Computer Modelling of Surfaces Fig. 6 Calcite (lOi0)thin film (25 A) with hydrogenphosphate ions (emphasised by the tetrahedra) (a) before and (b) after relaxation showing the formation of voids in the surface region preventing the crystal aggregating.This disruption of the crystal surface caused by the hydrogenphosphate additives shows the potency of the phosphate ions on the nucleation of calcite crystallites as observed experimentally. Other factors which can have a considerable impact on morphology include the degree of supersaturation, and variation in pH and ionic strength. In recent work2' Hopwood and Mann at Bath have systematically varied each of the properties in the study of the growth and morphology of barium sulfate crystals. Importance of Surface Charge to Barium Sulfate The importance of ionic strength in crystal growth is that when increased significantly it will stabilise charged surfaces. This occurs because on increasing ionic strength the charge compensatory space charge layer increases in density as it becomes thinner.This will have the effect of increasing the cohesion and hence the surface stability. There will also be a secondary effect of allowing easier access of like ions of the same charge as the surface and in the extreme case may allow crystallites to join. We have performed some preliminary calculations on BaSO,, and have attempted to correlate the simulation results with the morphological changes observed by Hopwood. We leave the detailed description of the structure and surfaces of barium sulfate to Allan and Rohl (this Discussion). We found that using the less sophisticated potential the most stable surfaces were still predicted to be the (210)and (001) surfaces; the energies are given in Table 4.However, to estimate the effect of high ionic strength on morphology we need to estimate which surfaces can support a high surface charge. Frenke122 first developed an approach for calculating surface charge from the S. C. Parker et al. differences in free energies of formation of surface anion and cation vacancies. Recently this was extended to other ceramic materials by Duffy and Ta~ker,~, who also showed that the surface charge can be calculated from the difference in vacancy interaction energies (DVIE). This is simply the difference in the segregation energy of a barium vacancy compared to the segregation energy of a sulfate vacancy. These differences, neglecting relaxation around the defects, are given in Table 4.The larger the magnitude of DVIE the larger is the surface charge; a negative value implies a negatively charged surface. Neglecting the (101) and { 11 l} surfaces the most highly charged surfaces are predicted to be (011) followed by (loo}. These two surfaces are observed to dominate the morphology at high ionic strength. The effectiveness of the space charge depends on the valence of the ions contributing to it, the charged surfaces are stabilised by adding more highly charged ions. For example, the presence of multiply charged anions should further stabilise the positively charged surfaces, namely the (01 1) and the (100) surfaces. The reason for neglecting the { lOl} and (1 1 1) is based purely on kinetic arguments.In each case there are two possible terminations with similar surface energies (denoted in Table 4 by the suffix ‘a’ and ‘b’ after the index). On the macroscopic scale we would expect these surfaces to be composed of regions of both terminations, and hence they will be atomistically ‘rough’. This roughening of the surface should lead to rapid growth since step sites are centres for growth. The correlation implies that further work including the full defect relaxation and solving the space charge problem explicitly may allow us to make quantitative predictions. Conclusion We have reported current work on the prediction of the effect of surfaces defects and additives on the morphologies of NiO and calcite using atomistic simulation.These predictions are consistent with experimental observations and demonstrate the potential of atomistic simulations in modelling the surface structure and growth of inorganic materials under a range of conditions. As with all simulations, the reliability and hence the confidence we can place in the results will continue to improve as we develop more sophisticated interatomic potentials. For example, improvement in the derivation of interatomic potentials for CaCO, will enable us to extend the range of applications. This includes determining the factors which control calcite formation over that of the other polymorphs, i.e. aragonite and vaterite. Further developments currently under investigation include the effect of temperature, particularly for ceramics such as NO, where temperatures can exceed 1500 K.In addition, studies are currently in progress to model the effect of pH, ionic strength and solvated ions on the surface stability. The preliminary results on the effect of phosphate on calcite nucleation and the correlation between surface vacancy formation energies of free surfaces and ionic strength suggest that atomistic simulation can aid in the interpretation, and perhaps the prediction, of experimental growth data. Table 4 Surface energies and differences in vacancy interaction energies (DVIE) surface surface energy/J m-’ DVIE/eV 210 0.48 1.o 001 0.61 -1.0 010 0.69 -1.2 llla 0.76 3.4 lllb 0.86 2.2 lOla 0.78 -6.1 lOlb 0.95 5.0 100 0.84 1.4 01 1 0.98 5.0 84 Computer Modelling of Surfaces References 1 P. W.Tasker, E. A. Colbourn and W. C. Mackrodt, J. Am. Ceram. Soc., 1985,68, 74. 2 P. W. Tasker and D. M. Duffy, Surf Sci., 1984, 137,91. 3 N. L. Allan, P. Kenway, W. C. Mackrodt and S. C. Parker, J. Phys. Condens. Matter, 1989, 1, SB119. 4 P. J. Lawrence and S. C. Parker, in Computer Modelling of Fluids, PolymersandSolids, ed. C. R. A. Catlow, S. C. Parker and M. P. Allen, Kluwer, Dordrecht, 1990, vol. 293, pp. 219-248. 5 M. J. L. Sangster and A. M. Stoneham, Philos. Mag., 1981,43, 597. 6 R. A. Jackson and G. D. Price, Mol. Simul., in the press. 7 F. Bertaut, Compt. Rend., 1958,246, 3447. 8 G. Wulff, Z. Kristallogr. Kristallgeom., 1901, 34, 949. 9 J. W. Gibbs, Collected Works, Longman, New York, 1928.10 P. M. Oliver, S. C. Parker and W. C. Mackrodt, Modelling Simul. Muter. Sci. Eng., submitted. 11 C. A. Handwerker, M. D. Vaudin and J. E. Blendell, J. Phys. (Paris) C.5, 1988,49, 367. 12 P. W. Tasker and D. M. Duffy, Philos. Mag. A, 1986,54, 759. 13 M. P. C. Weijnen and G. M. van Rosmalen, Industrial Crystallisation, 1984, 84, 61. 14 S. Mann, S. J. Didymus, N. P. Sanderson, N. P. Heywood and E. J. Aso-Samper, J. Chem. Soc.,Faraday Trans., 1990, 86, 1873. 15 H. Sawada and Y. Takuechi, 2. Krystallogr., 1987, 181, 179. 16 J. M. Didymus, Ph.D Thesis, University of Bath, 1992. 17 J. 0.Titiloye, S. C. Parker, D. J. Osguthorpeand S. Mann, J. Chem. Soc., Chem. Commun., 1991,20,1494. 18 H. Leitmeier, Neues. Jahrb. Miner. Abh., 1910, 1, 49. 19 S. Rajam and S. Mann, J. Chem. Soc., Chem. Commun., 1990, 1789. 20 T. Suzuki, S. Inomata and K. Sawada, J. Chem. Soc., Faraday Trans., 1986,82, 1733. 21 J. D. Hopwood, J. Crystal Growth, submitted. 22 J. Frenkel, in Kinetic Theory ofLiquids, OUP, New York, 1946, p. 37. 23 P. W. Tasker and D. M. Duffy, Philos. Mag. A, 1984, 50, 143. 24 R. A. Jackson and P. Meenan, in the press. 25 B. W. H. van Beest, G. J. Kramer and R. A. van Santen, Phys. Rev. Lett., 1990,64, 1955. 26 P. Saul, C. R. A. Catlow and J. Kendrick, Philos. Mag. B, 1985, 51, 107. Paper 3/00227F; Received 1 1 th January, 1993