Molecular dynamics simulations of Na deposition on the TiO2(110) surface

Surface Science 409 (1998) 92–100

Molecular dynamics simulations of Na deposition on the TiO (110) surface 2 M.A. San Miguel, C.J. Calzado, Javier Fdez. Sanz * Departamento de Quı´mica-Fı´sica, Facultad de Quı´mica, Universidad de Sevilla, E-41012 Seville, Spain Received 11 November 1997; accepted for publication 17 March 1998

Abstract A computational study of the sodium deposition on the TiO (110) rutile surface was carried out. Using both molecular dynamics 2 methods and complementary ab initio Hartree–Fock embedded cluster calculations, the Na arrangements on the surface for coverages up to h=0.5 ML were considered. Our simulations show that Na adsorbs on the surface binding three oxygen atoms, two protruded and one basal, in a quasi-symmetric configuration. This coordination type is found whatever the coverage, in agreement with previous experimental and theoretical works. The Na atoms are found to be arranged in an alternated disposition around the protruded oxygen atom rows. This distribution is almost perfect for a h=0.25 ML coverage, since coadsorption effects are small and the Na atoms find place enough to be accommodated. For a coverage of h=0.5 ML, some islands, in which the Na atoms are facing each other across two protruded oxygen atoms, are found, even in long time simulations. Formation of Na O dimers was not 2 observed. © 1998 Elsevier Science B.V. All rights reserved. Keywords: Ab initio quantum chemical methods and calculations; Adatoms; Alkali metals; Models of surface chemical reactions; Molecular dynamics; Surface structure; Titanium oxide

1. Introduction Metal deposition on a variety of substrates constitutes one of the most useful methods to obtain modified materials which exhibit appealing technological properties. Among a wide variety of supports, metal oxides are often employed because of their temperature and mechanical resistance. One of the major applications of metal deposition on metal oxides refers to the promotion effects within heterogeneous catalysis. For this purpose, alkali metals are of special interest since, because of their low first ionization potentials, they behave as efficient surface reductors. In spite of the impor* Corresponding author. E-mail: [email protected] 0039-6028/98/$19.00 © 1998 Elsevier Science B.V. All rights reserved. PII: S0 0 39 - 6 0 28 ( 98 ) 0 02 4 5 -3

tance of these phases, the detailed structural and electronic modifications induced on both the substrate and the metal after deposition are not wellknown, and the literature concerning these systems is relatively scarce. Because of both simplicity and catalytic interest, rutile TiO promotion by deposition of Na atoms 2 constitutes one of the most studied phenomena [1–3]. As shown by Onishi et al. [4], Na deposition onto the (110) face of rutile TiO induces a 2 lowering of the work function of 3.4 eV. The appearance of emission in the bulk bandgap in UPS spectra indicated that Na was ionically adsorbed while the Ti4+ surface cations were reduced to Ti3+. From the observed c(4×2) LEED pattern at 0.5 ML Na coverage, these

93

M.A. San Miguel et al. / Surface Science 409 (1998) 92–100

authors suggested that Na adsorbs, forming Na O dimers in which the Na+ ion is coordinated 2 to two protruded surface oxygen atoms and one basal oxygen atom in a symmetric arrangement. In a further work, Murray et al. [5] reported scanning tunneling microscopy (STM ) experiments conducted on reduced (1×2) TiO (110) 2 samples and low Na coverage. Based on the STM images, these authors proposed a new model in which the Na atoms were adsorbed on an adjacent position to only one bridge oxygen. However, in a more recent work, Nerlov et al. [6 ] showed that the two different adsorption models correspond in fact to the different structures of the surfaces, and that the (1×2) reconstruction of the reduced (110) TiO would favor the adsorption adjacent 2 to the protruded oxygen, while in the stoichiometric surface, the adsorption between two bridge oxygen atoms would be preferred. The electronic details concerning the interaction between atomic Na and the TiO rutile surfaces 2 have been the subject of recent theoretical work in our group [7,8]. As shown by our ab initio embedded-cluster HF and CASSCF calculations, when a Na atom approaches the surface, an electron transfer from the 3s Na orbital towards the 3d Ti orbitals occurs, giving rise to Na+ species ionically adsorbed on the surface. As far as the (110) surface is concerned, calculations show that such a reduction takes place selectively on the penta-coordinated Ti centers, in agreement with the data reported by Onishi et al. [4]. Moreover, optimization of the Na position on the surface leads to an equilibrium geometry in which the Na+ ion is bound in a symmetric way to two protruded bridge oxygens and to the basal oxygen. This structure was found to agree quite satisfactorily with the model proposed by Onishi et al. [4], and the interpretation reported by Nerlov et al. [6 ]. However, because of the limited size of the clusters, the effects arising from the interaction between adsorbed Na atoms could not be included. Nevertheless, the size of these interactions is expected to be large, mainly for high coverages, and the question of whether such coadsorption effects could modify the description still remains open. In the present work, a molecular dynamics (MD)

simulation of the Na deposition onto the rutile (110) TiO surface is reported. Since our goal is 2 to analyze the dependence of the process on the surface Na–Na interactions, these MD simulations have been carried out for different Na coverages. On the other hand, these simulations have been complemented with ab initio embedded-cluster Hartree–Fock (HF ) calculations on selected structures with the aim of supporting some aspects of the MD simulations, in particular those dealing with the formation or not of surface Na O dimers. 2 The paper is arranged as follows. In Section 2, the models used to represent the rutile (110) surface, as well as a summary of the computational methods concerning both MD simulations and HF calculations, are reported. Section 3 is devoted to the results and discussion, while in Section 4 the main conclusions are outlined.

2. Surface models and computational methods To simulate TiO rutile, a system consisting of 2 1620 particles (540 Ti+1080 O) was arranged into a tetragonal computational box with a=26.63, b= ˚ . Free surfaces of TiO were 32.485, c=19.491 A 2 simulated through a slab built up by imposing periodic boundary conditions on this box in the three directions, but leaving a vacuum space of ˚ between the slabs in the direction perpendic100 A ular to the (110) plane. Addition of Na atoms was performed on both sides of the slab. The interaction between particles was described by pair potential functions consisting, as usual, of Coulombic and short-range terms: qq V = i j +V (s−r) ij ij r ij where q and q are the effective charges and r i j ij the interparticle distance. For the interaction between Ti and O atoms, a Buckingham type short-range potential was used:

A

B

A +A −r CC j ij V (s−r)=− i j +f (B B ) exp i i j ij r6 B +B ij i j where the parameters A, B, C and f were those given by Matsui and Akaogi [9]. For the inter-

94

M.A. San Miguel et al. / Surface Science 409 (1998) 92–100

action between Na+ and surface Ti and O atoms, a Pauling type potential was chosen:

A

B

s +s n 1 i j V (s−r)= ij n(s +s ) r i j ij where s , s are the ionic radii for species i and j i j [10] and the parameter n was set to 9, according to Adams and McDonald [11]. The whole set of charges and ionic radii for this potential are summarized in Table 1. According to Matsui and Akaogi [9], the effective charges of Ti and O atoms for unreduced TiO were set to +2.196 e and −1.098 e. Once the 2 surfaces have been reduced, the Ti charge has to diminish; however, as shown both experimentally and theoretically, such a reduction only takes place on the penta-coordinated surface Ti centers and involves an almost complete charge transfer of one electron from the Na atoms. Therefore, it seems reasonable to assume that upon reduction, the new species on the surface are Na+ cations (q=1 e) and the same number of reduced penta-coordinated Ti atoms to which an effective charge of +1.196 e can be ascribed. This assumption allows us to preserve the short-range contribution to the Ti–O pair potential, and to distinguish the Ti(III ) centers only through the coulombic term. MD simulations were performed in the microcanonical ensemble using the SIMULA code [12]. Numerical integration was carried out using the leap-frog algorithm with a time step of 1 fs. The temperature of the MD simulations was 300 K. The surfaces were thermally equilibrated in a run of 10 ps, in which the velocities were rescaled. After this thermalization period, another run of 5 ps was performed in order to assure that the system reached stability at 300 K. Finally, a proTable 1 Summary of the Pauling pair potential parameters used in the MD simulations

Ti(III ) Ti(IV ) O Na

q (e)

˚) s (A

1.196 2.196 −1.098 1.0

0.68 0.68 1.40 0.95

duction run of 10 ps was performed in order to calculate the statistics of the system. In order to analyze the interaction between sodium atoms and the surface using a finite approach, a set of ab initio Hartree–Fock embedded cluster calculations was performed. The cluster models to describe the surface were carefully described in Ref. [8] and only a brief summary will be reported here. Basically, the surface was described by a cluster of formula Ti O16− sur4 16 rounded by a set of 17 total ion potentials representing the nearest Ti ions directly bound to the cluster oxygen atoms. The whole system was then embedded in an array of 1404 point charges to introduce the Madelung potential of the crystal at the surface. In this cluster model, both five-fold and six-fold coordinated Ti atoms are present. The core electrons of Na and Ti atoms (10 electrons) were described by using the effective core potentials proposed by Hay and Wadt [13– 15]. The basis sets for the valence electron were (10s5p5d )/[2s2p2d ] for Ti and (3s1p)/[2s1p] for Na. For oxygen atoms, the (9s5p)/[4s3p] basis set optimized by Broughton and Bagus was used [16,17]. For the TIP representation of Ti ions, the large core effective potentials of Hay and Wadt [13–15] were chosen. Optimized geometries were obtained from unrestricted HF wave functions using standard analytical techniques for gradient computation. All the calculations were carried out using the HONDO8.4 program [18] running on a DEC ALPHA 4100 computer.

3. Results The first step in the study was to simulate the isolated (110) surface. After equilibration, the surface was found to relax involving displacements of the first few outlayers. In general, the surface ˚ for the atoms were found to move inward: 0.1 A ˚ for the pentaprotruded oxygen ions, and 0.26 A coordinated Ti atoms. These values are in agreement with previously reported results obtained from periodic HF (Reindhart and Hess [19]), periodic LDA (Ramamoorthy et al. [20]), and molecular mechanics (Oliver et al. [21]), as

95

M.A. San Miguel et al. / Surface Science 409 (1998) 92–100

Table 2 Summary of surface relaxation obtained from MD simulations and comparison with other published data. Displacements of surface ˚ ions with respect to their ideal crystal positions are in A Iona

Ti p

Ti h

O br

O ba

O h

O p

MD (this work) Periodic HFb Periodic LDAc Mol. mech.d Exp.e

−0.26 −0.15 −0.17 −0.25 −0.16±0.05

0.02 −0.09 0.13 0.25 0.12±0.05

−0.11 −0.14 −0.07 −0.08 −0.27±0.08

0.0 0.07 0.13 0.02 0.05±0.05

0.0 −0.07 −0.08 — 0.05±0.08

0.0 −0.02 0.02 — 0.0±0.08

a Ti and Ti refer to penta- and hexa-coordinated surface Ti. O =bridging protruded oxygens, O =basal (second oxygen layer). p h br ba O =oxygen bound to Ti (third oxygen layer). O =oxygen bound to Ti (fourth oxygen layer). h h p p b Ref. [19]. c Ref. [20]. d Ref. [21]. e Ref. [22].

well as with experimental data recently reported by Charlton et al. [22]. For the sake of comparison, these data are collected in Table 2, and a general view of the surface is reported in Fig. 1. The Na deposition was then simulated by addition of a given number of Na+ cations to the surface while, according to the above discussion, the same number of penta-coordinated Ti atoms were formally reduced from Ti(IV ) to Ti(III ). Four series of simulations were performed, corre-

Fig. 1. General view of the TiO (110) rutile surface as obtained 2 from molecular dynamics simulations.

sponding to the addition of (i) one Na atom, (ii) two Na atoms, (iii) 0.25 ML and (iv) 0.5 ML. The addition of only one or two Na atoms to the surface by computational box can be considered as unrealistic since it corresponds to really low coverages; however, these systems are of some interest in the present context since they allow for a meaningful comparison with the results obtained from finite ab initio HF calculations. On the other hand, the definition of coverages follows that of Onishi et al. according to which a coverage of 0.5 ML corresponds to the deposition of as many Na atoms as penta-coordinated Ti atoms. Results for adsorption of a single atom are reported in Fig. 2 and Table 3. As shown in the snapshot of the figure, the Na+ ion is found to bind two bridging oxygens and one basal oxygen atom reaching a three-fold coordination. The Na+ ion lies somewhat above the surface, and in fact it appears to swing around the bridge O–O axis, although the averaged position observed is that in which the Na ion lies by the side of the Ti(III ) center, according to the lower Ti(III )–Na coulombic repulsion. This structure agrees with the models proposed from experiment by Onishi et al. [4], and Nerlov et al. [6 ], as well as with that found from ab initio HF calculations [8]. As can be seen in Table 3, HF and MD geometrical parameters are found in excellent agreement, the larger deviation being the distance between Na ˚ shorter in and the basal oxygen, Na–O , 0.14 A ba the MD simulations. In the same table, and for

96

M.A. San Miguel et al. / Surface Science 409 (1998) 92–100

( Fig. 3, top). In the second possibility, the two Ti ions belong to different rows, i.e. they are p separated by a row of protruded oxygen atoms ( Fig. 3, bottom). For the first possibility, the most suitable disposition appears to be that in which the Na+ ions are alternated around the same Ti p row, but interacting with oxygen ions of different rows, as shown in the top part of Fig. 3. For the second possibility, one always finds a Na+ ion near a Ti , while the other one may occupy several p positions, one of them is shown in the bottom part of Fig. 3. In any case, it is worth noting that the coordination sites are always of the same type, with Na+ ions binding quasi-symmetrically two

Fig. 2. MD simulations snapshots for the deposition of a single Na atom onto the TiO (110) rutile surface. 2 Table 3 Structural parameters for the adsorption of a single Na atom on the TiO (110) surface obtained from ab initio embedded2 cluster Hartree–Fock calculations and MD simulations. Bond ˚ distances in A

d(Na–Ti ) p d(Na–Ti ) h d(Na–O ) br d(Na–O ) ba

HFa

MD

b

4.54 3.05 2.36 3.38

4.50 3.07 2.32 3.24

3.23 2.68 2.29 2.29

a Ref. [8]. b Ref. [6 ].

the sake of comparison, the distances proposed by Nerlov et al. [6 ] from a simple geometrical model and the ionic radii are also reported. When two Na atoms are added to the surface, two main cases can be considered. In the first one, the reduction involves non-contiguous five-fold coordinated Ti atoms, Ti , in which case, the p interaction between the associated Na+ ions is small, and the situation is similar to that found in the previous paragraph. In the case that the two Ti are close each to the other, two possibilities p can be distinguished. In the first one, they are contiguous, i.e. they fall in the same row of Ti p

Fig. 3. MD simulations snapshots for the deposition of two Na atoms onto the TiO (110) rutile surface. 2

M.A. San Miguel et al. / Surface Science 409 (1998) 92–100

bridging oxygens and, in less extension, a basal one. Results for larger coverages are reported in Fig. 4. The top part of the figure corresponds to h=0.25 ML, and the bottom to h=0.5 ML. As can be seen, the Na surface coordination scheme drawn above remains almost unchanged for h= 0.25 ML, with Na+ ions binding symmetrically two bridging oxygens and a basal one. Also, a general trend for the Na+ ions to alternate on each side of an oxygen atom row is observed. For a coverage of h=0.5 ML, most of the Na+ ions are also found with the same type of coordination, and the formation of alternated Na chains appears

Fig. 4. MD simulations snapshots for the deposition of Na onto the TiO (110) rutile surface. Top, h=0.25 ML; bottom, h= 2 0.5 ML. The cage highlights the less stable arrangement.

97

sharply. However, the presence of some islands with a different Na arrangement type is clearly observed. In these regions, the Na+ ions appear to be facing each other across a pair of oxygen atoms of the same row. From a purely electrostatic point of view, such an arrangement is less stable since it involves a larger repulsion between the Na+ ions, however, it appears to persist even for larger simulation times. Another interesting result of our simulations refers to the absence of Na+ ions adsorbed in an adjacent arrangement. In this disposition, the Na+ ions would be assembled in pairs, one on each side of the same protruded oxygen atom. This result confirms the interpretation of Nerlov et al., in the sense that they appear only in the modified (110) surface. In order to obtain an estimate of the relative stability for the three different arrangements, a set of ab initio embedded cluster HF calculations for the cluster described in the computational section and two Na atoms was carried out. The three models, labeled a, b and c, are depicted in Scheme 1 (see also Figs. 5 and 6). Model a corresponds to the alternated disposition of Na+ cations along the bridging protruded oxygen row. Model b is that called adjacent, while in model c the Na+ ions are facing between two protruded oxygen ions. In these calculations, the position of the two Na+ ions was optimized assuming no relaxation of the surface. The results are reported in Table 4 and Figs. 5 and 6. As can be seen in the table, the most stable structure corresponds to model a, about 25 kcal/mol lower than the adjacent position b. Model a is a true minimum in the potential energy hypersurface, while model b is unstable with respect to model a, and can only be obtained by imposing C symmetry. On the 2v other hand, we were unable to locate a stationary point for model c, since it always evolves towards

Scheme 1. Arrangements of two Na atoms round a row of protruded oxygens.

98

M.A. San Miguel et al. / Surface Science 409 (1998) 92–100

Fig. 5. Optimized structure for model a obtained from ab initio HF embedded cluster calculations. Top, view perpendicular to the (110) plane.

model a. Comparing the structural parameters with those of Table 3, it turns out that Na+ cations are now closer to the surface, suggesting that some coadsorption effects are in play. Thus, the distance between Na and five-fold coordinated Ti ions p ˚ , while the distance decreases from 4.54 to 3.88 A between Na and hexa-coordinated Ti centers h˚ remains almost unchanged (3.05 vs. 2.94 A ). These displacements, which in fact almost correspond to a rotation around the Ti centers, are due to the h repulsion between the Na+ ions across the oxygen rows. In Table 4, the values for these parameters obtained from the MD simulations are also reported. These values were obtained from the peaks of the radial distribution functions g(r) of Fig. 7. As shown, the quoted interatomic distances

Fig. 6. Optimized structure for model b obtained from ab initio HF embedded cluster calculations. Top, view perpendicular to the (110) plane. Table 4 Structural parameters for the adsorption of Na on the TiO (110) surface obtained from ab initio embedded-cluster 2 Hartree–Fock calculations and MD simulations. Bond dis˚ and relative energies in kcal/mol tances in A HF

d(Na–Ti ) p d(Na–Ti ) h d(Na–O ) br d(Na–O ) ba d(Na–Na) E rel

MD

a

b

a

c

3.88 2.94 2.38 2.79 3.76 0

3.15 3.46 2.23 2.94 3.60 25

3.90 2.96 2.40 2.45 3.90 —

3.50 3.01 2.40 2.41 2.40 —

M.A. San Miguel et al. / Surface Science 409 (1998) 92–100

99

However, once such irregular domains are present, the system can hardly reach a perfect alternated arrangement since it would involve a concerted displacement of several rows of Na+ ions. That is why such domains persist even for long time MD simulations.

4. Conclusions

Fig. 7. Partial radial distribution functions (arbitrary intensities) obtained from MD simulations. Top, h=0.25 ML; bottom, h=0.5 ML.

are in rather good agreement with those predicted by the HF calculations. There is only a noticeable exception: the distance between Na and the basal oxygen, d(Na–O ), which is found to be signifiba ˚ cantly shorter (0.34 A ) in the MD calculations. This lowering is mainly due to two reasons: the first is the surface relaxation induced by Na deposition, involving a displacement outwards of the surface of the basal oxygen ion; the second is the larger separation between Na+ ions across the oxygen rows because of their repulsion. The final situation is that in which the Na–O interatomic distances are almost the same. The relative stability of these models allow us to confirm that the alternated arrangement mostly observed in the MD simulation snapshots of Fig. 4 results simply from a minimization of the Na+–Na+ repulsion, and that the small islands observed correspond to more energetic situations.

In this work we report on a molecular dynamics simulation of the Na deposition at the TiO (110) 2 rutile surface. According to both experimental and theoretical results, the simulations have been carried out under the assumption that, for moderate coverages, deposition of Na involves a selective reduction of five-fold coordinated Ti ions. Simulations at very low coverages unambiguously show that the preferred coordination site corresponds to an adsorption ‘‘in between’’ the bridging (protruded ) oxygen atoms, in which the Na+ ion binds two protruded oxygen atoms and a basal one in a quasi-symmetric configuration. On these sites, the Na+ ions lie somewhat above the surface, but closer to the reduced Ti center. These results agree satisfactorily with previous theoretical calculations. For higher coverages, the coordination ambient is found to be essentially the same, but with the Na atoms closer to the basal oxygens. Also, the Na atoms are found to be arranged in an alternated disposition around the protruded oxygen rows. For the larger coverage analyzed here, h=0.5 ML, which corresponds to as many Na atoms as fivefold coordinated Ti atoms, some domains with a different Na arrangement are observed. In these, the Na+ ions have the same coordination but they are faced across a pair of protruded oxygen atoms. On the other hand, adsorption in sites adjacent to the protruded oxygen atoms is not observed, confirming the interpretation of Nerlov et al. [6 ]. Finally, in disagreement with the proposal of Onishi et al. [4], formation of Na O dimers is not 2 found at these coverages. Ab initio Hartree–Fock embedded cluster calculations carried out on some representative models show that the alternated arrangement is by far the most stable. Moreover, the structural parameters

100

M.A. San Miguel et al. / Surface Science 409 (1998) 92–100

optimized at this level of calculation are found in quite good agreement with the data obtained from the radial distribution functions arising from MD simulations. These results support the reliability of the pair potentials used here. From both the MD simulations and the HF calculations, the Na deposition on this surface can be described as the Na+ ions mostly occupying three-coordinated sites following an alternated arrangement around the oxygen protruded rows. However, such an arrangement is not perfect since the random nature of the deposition leads to formation of local irregularities. Also, some irregular domains in which the Na+ ions are facing each other across a pair of oxygen atoms of the same row may appear. Although both types of defect correspond to more energetic arrangements, they can hardly reach an ideal alternated configuration since it would involve a simultaneous displacement of entire Na rows.

Acknowledgements This work was supported by the DGICYT (Spain, Project No. PB95-1247), and by the European Commission (Contract No. ERBCT1CT94-0064). Part of this work was presented at the 9th International Congress of Quantum Chemistry held in Atlanta, 9–14 June 1997. The authors are indebted to the participants of this meeting for many helpful suggestions.

References [1] V.E. Henrich, P.A. Cox, The Surface Science of Metal Oxides, Cambridge University Press, Cambridge, 1994.

[2] P.A. Cox, Transition Metal Oxides, Clarendon Press, Oxford, 1995. [3] K. Prabhakaran, D. Purdie, R. Casanova, C.A. Muryn, P.J. Hardman, P.L. Wincott, G. Thornton, Phys. Rev. B 45 (1992) 6969. [4] H. Onishi, T. Aruga, C. Egawa, Y. Iwasawa, Surf. Sci. 199 (1988) 54. [5] P.W. Murray, N.G. Condon, G. Thornton, Surf. Sci. 323 (1995) L281. [6 ] J. Nerlov, S.V. Christensen, S. Weichel, E.H. Pedersen, P.J. Mo¨ller, Surf. Sci. 371 (1997) 321. [7] C.J. Calzado, J. Oviedo, M.A. San Miguel, J.F. Sanz, J. Mol. Catal. 119 (1997) 135. [8] M.A. San Miguel, C.J. Calzado, J.F. Sanz, Int. J. Quantum Chem., in press. [9] M. Matsui, M. Akaogi, Mol. Sim. 6 (1991) 239. [10] L. Pauling, The Nature of the Chemical Bond, Cornell University Press, Ithaca, NY, 3rd edn., 1960, p. 514. [11] D.J. Adams, I.R. McDonald, Physica B 79 (1970) 159. [12] SIMULA – A molecular dynamics and visualization software developed by L.J. Alvarez, M.A. San Miguel, Universidad Nacional Auto´noma de Me´xico, 1995. [13] P.J. Hay, W.R. Wadt, J. Chem. Phys. 82 (1985) 270. [14] P.J. Hay, W.R. Wadt, J. Chem. Phys. 82 (1985) 284. [15] P.J. Hay, W.R. Wadt, J. Chem. Phys. 82 (1985) 299. [16 ] J.Q. Broughton, P.S. Bagus, Phys. Rev. B 30 (1984) 4761. [17] J.Q. Broughton, P.S. Bagus, Phys. Rev. B 36 (1987) 2813. [18] M. Dupuis, S. Chin, A. Ma´rquez, CHEM-Station and HONDO: modern tools for electronic structure studies including electron correlation, in: G.L. Malli ( Ed.), Relativistic and Electron Correlation Effects in Molecules and Clusters, NATO ASI Series, Plenum Press, New York, 1994. [19] P.B. Reinhardt, A. Hess, Phys. Rev. B 50 (1994) 12015. [20] M. Ramamoorthy, D. Vanderbilt, R.D. King-Smith, Phys. Rev. B 49 (1994) 16721. [21] P.M. Oliver, G.W. Watson, E.T. Kelsey, S.C. Parker, J. Mater. Chem. 7 (1997) 563. [22] G. Charlton, P.B. Howes, C.L. Nicklin, P. Steadman, J.S.G. Taylor, C.A. Muryn, S.P. Harte, J. Mercer, R. McGrath, D. Norman, T.S. Turner, G. Thornton, Phys. Rev. Lett. 78 (1997) 495.