Conference-ICSFS-14-Testing Interatomic Potentials for QM / MM Embedded-Cluster Calculations on Ceria Surfaces ∗

Interatomic potentials (IP) have demonstrated considerable application in the study of bulk and surfaces of ceria (CeO2), and also the long range component in QM/MM calculations. Despite the development and ability of several IPs to reproduce the structural and dielectric properties of bulk ceria, in the absence of high quality electronic structure calculations it remains unclear how well these models perform when applied to surfaces. In this paper we present a comparison of several IPs from the literature in comparison with periodic density functional theory (GGA+U) for the calculation of the (111), (110) and (100) low index surfaces of ceria. While the IP approach reproduces the correct order of surface stability ((111)> (110) > (100)) and relaxation ((100)> (110) > (111)) found in the application of GGA+U , the surface energies calculated with IPs are considerably larger than those found with GGA+U . Based on the most comparable IP tested against GGA+U , we have performed trial QM/MM embedded-cluster calculations on the stoichiometric and defect (110) surface of ceria and compared the resulting defect energy with GGA+U . [DOI: 10.1380/ejssnt.2009.413]


I. INTRODUCTION
Ceria (CeO 2 ) is an insulating, non-magnetic rare-earth oxide [1]; consisting of a cubic fluorite structure with four cerium and eight oxygen atoms per unit cell and an experimental bulk lattice parameter of 5.411 Å [2].The large number of experimental and theoretical studies of ceria in the literature attests to its technological importance, especially its role in automobile three-way catalytic converters.The reversible formation of partially reduced CeO 2 surfaces, creating an oxygen vacancy and two neighbouring formally charged Ce 3+ ions, occurs due to the removal of an oxygen atom from a surface lattice site.The reduction process is more favourable on ceria, compared to other oxides, allowing facile oxidation/reduction of environmentally harmful molecules such as CO and NO 2 at the temperatures present in an exhaust gas [3].
A number of experimental studies concerning low index stoichiometric [3][4][5][6][7] and reduced [8] ceria surfaces exist in the literature.From scanning tunneling microscopy (STM) [3], ion scattering spectroscopy [6] and atomic force microscopy (AFM) [7] techniques it can be observed that the (111) surface is oxygen terminated with small ionic relaxations [3].The (110) surface is terminated with a stoichiometric CeO 2 anion-cation-anion atomic plane, as described by STM and electron diffraction methods [4,5].The least stable surface is the (100) surface, and upon cleaving a dipole moment is produced perpendicular to the surface.Strong surface reconstruction is required to cancel out the respective dipole moment through the removal of 50% of surface terminating oxygen atoms [5,6].Of the possible terminations of the (100) surface, the oxy-gen terminated surface has the lowest energy [5,6], and is accepted to be the experimentally observed termination [5].Studies of partially reduced ceria surfaces have been carried out using several experimental techniques and are described in a recent review concerning oxygen vacancies in metal oxides by Ganduglia-Pirovano et al. [8].The main observations from experimental research on reduced ceria surfaces conclude that Ce 3+ ions can be found in surface and sub-surface layers, and that the +4 to +3 reduction of ceria is driven by the excess electrons that remain (occupying strongly localized f states) after oxygen removal occurs.
The aim of our current research is to study the re- duction of ceria surfaces (oxygen vacancy defects) using the hybrid quantum mechanics / molecular mechanics (QM/MM) embedded-cluster approach presented by Sokol et al. [31].However, it is not clear which IP gives rise to the most accurate model of ceria.Although all of the available IP models can reproduce experimental data for bulk ceria (since they are fitted to data for the bulk), it is not clear how accurate these models are for the study of surface structures and energetics.One method to asses the applicability of IPs available for ceria in the literature is to use ab initio calculations as a benchmark with which to compare their suitability for surface calculations.While several publications exist that study individual IPs against ab initio ceria calculations (HF [23][24][25] and DFT [25,32]), no such wide ranging study as described above has been performed to date.Thus, as a pre-requisite to QM/MM calculations, we present a comparison of properties of the ( 111), ( 110) and (100) surfaces of ceria using a number of self consistently derived interatomic potentials against ab initio GGA+U .The choice of GGA+U is motivated in part by the body of data showing that it provides good results for the structure and energetics of many oxides [33] as well as for ceria [24,25,27,28].Following this, we present a QM/MM calculation on the reduction of the (110) ceria surface, utilizing one of the IPs tested to illustrate the applicability of our approach.

II. COMPUTATIONAL METHODS
The structures of the low index (111), (110) and (100) ceria surfaces are displayed in Fig. 1.The Tasker classification of these surfaces into type 1 (110), type 2 (111) and type 3 (100) has been discussed previously [34].For the bulk terminated type 3 (100) surface with a dipole moment perpendicular to the surface, the dipole is removed by transferring half of the surface layer from one side of the slab to the other, consistent with experimentally observed oxygen termination of the (100) surface previously described in the introduction, producing a 50% vacant oxygen surface layer [5,6].

A. Slab Surface Calculations
Atomistic simulations were carried out with the molecular modeling (MM) package GULP [35], in which the total lattice energy of a periodic system is computed as a sum of the energies of long-range and short-range inter-actions, where r ij is the interatomic distance, q i and q j are the ionic charges and A ij , ρ ij and C ij are parameters for the short range interactions, with a 20 Å cut-off.Ionic polarisability is treated by the shell model, which couples a massless shell of charge Y to a point core with a potential of the form, where (δr i ) is the distance between the shell and the core and k 1 and k 2 are the harmonic and quartic spring constants, respectively.The ionic charges for the potential models are taken from formal oxidation states, so that q Ce is +4.0 and q O is −2.0.All interatomic potential parameters used in this paper are presented in Table I.
The GGA+U calculations were performed using the Vienna ab initio simulation package VASP [36][37][38] which utilizes a plane-wave basis set approach.The projectoraugmented wave (PAW) [39,40] method was used to treat the valence core interaction with cores of Ce:[Xe] and O:[He].The generalized gradient approximation (GGA) was used in the form of the Perdew-Burke-Ernzerhof (PBE) [41] exchange-correlation functional.The U parameter in DFT+U was set to 5 eV [28,29] with a converged plane wave cutoff energy of 500 eV and a Brillouin zone grid sampling of 2×2×1 Monkhorst-Pack k-points.
Both GGA+U and IP surface calculations are performed using the same basic periodic slab, i.e. a crystal of particular finite thickness in a three dimensional periodic cell, which generates two surfaces via the introduction of a vacuum gap perpendicular to the surface.Demonstrating 4 terminating oxygen atoms, the (111) and (100) surfaces of ceria employ a p(2×2) expansion of the unit cell, while the (110) surface utilises a p(2×1) expansion.A 15 Å vacuum gap is large enough to remove interactions between periodic images perpendicular to the surface and the slab thickness is such that the structure in the middle of the slab is bulk-like.For the (111) surface, the slab thickness is 21/12 atomic layers (20.3/11.0Å), for IP/GGA+U calculations, respectively.The (110) surface contains 15/7 atomic layers (26.9/11.6Å) and the (100) surface consists of 29/9 atomic layers (37.7/10.4Å).Different slab thicknesses for the IPs and GGA+U calculations were utilized in order to maintain the minimum slab thickness required for surface energy convergence for each computational method.The surface energy is computed as E s = (E slab − E bulk )/2A, where A is the surface area and the factor of 2 arises from the presence of two surfaces in a slab model.Calculation of the respective surface energies allows a direct comparison of the IPs and GGA+U methods, with convergence to 0.01 J/m 2 in slab and vacuum thickness.

B. QMMM Calculations
The embedded-cluster quantum mechanics/molecular mechanics (QM/MM) approach is implemented in the computational chemistry software ChemShell [42], utilizing the approach described by Sokol et al. [31].Within ChemShell we have utilised GAMESS-UK [43] for the QM part of the calculations, while the MM section of the calculations used the GULP code [35].Our trial QM calculations on the (110) surface of ceria consisted of a Ce 7 O 30 /Ce 7 O 29 embedded-cluster and was calculated using the Hartree-Fock (HF) approach, along with effective core potentials (ECPs) for the Ce ions of the QM region [44] / QM interface region [45], and a Basis Set of DZV and DZVP for Ce [44] and O [46] ions, respectively.The hemi-spherical surface, constructed in the manner proposed by Sokol et al. [31], contains several enclosing spheres consisting of different regions of Ce and O ions, up to diameter of ∼50 Å.The central QM cluster of Ce and O ions (Ce 7 O 30 /Ce 7 O 29 ) is encased by 21 Ce QM interface ions, preventing the electrons from the QM region from collapsing into the MM region.These QM ions are surrounded by relaxed MM Ce and O ions up to a radius of 13 Å from the centre.The remaining sphere consists of Ce and O ions fixed in space using MM, and are surrounded by point charges to reproduce the electrostatic potential of the surface.The defect energy is calculated as , where E(CeO 2 ) is the total energy from the combined QM/MM calculation of the stoichiometric surface, E(CeO 2red ) is the total energy of the reduced surface and E(O 2gas ) is the QM calculated energy of an oxygen molecule.

A. Comparison of interatomic potentials
Three main groups of interatomic potentials (IPs) were tested and compared with density functional theory (GGA+U ) for the optimization of the low index (111), (110) and (100) surfaces of ceria and the calculation of their surface energies and geometries.The first group (interatomic potentials 1-3, IPs 1-3) contains a Ce core which is negatively charged and a shell that is positively charged; these shell model charges are "non-traditional".The first potential (IP1) is the Butler potential [13], in which the O 2− -O 2− parameters were derived by Catlow [47] and the Ce 4+ -O 2− interactions were derived using the electron gas method [48,49].IP2 is a modified version of IP1, where the O 2− -O 2− shell parameters were derived empirically by Lewis and Catlow [50], and the dispersion term (C ij ) for the Ce 4+ -O 2− interaction has been added along with modified k 1 terms and a new quartic k 2 term preventing the unphysical possibility that the shell will fall off (very large core-shell separation).IP3 is identical to IP2, with the exception of the lack of a k 2 term.Other theoretical research groups have also used identical or similar potentials to IPs 1-3 to study low index ceria surfaces [9][10][11][12][14][15][16].The second group, containing only the Grimes and Catlow [17] potential IP4, is fitted using the electron gas method of Wedepohl [48] for the A ij and ρ ij parameters of the Buckingham potential, while the dispersion term (C ij ) is calculated individually and incorporated directly into the potential.The third group (IP5 and IP6) [19] was derived to model high index surfaces such as (322), (311) and (331), based on an O 2− -O 2− potential whose accuracy in modeling surface and defect structures was previously demonstrated [51,52].The latter two groups have "traditional" shell model charges, with a positively charged core and negatively charged shell.

Bulk Properties
The lattice constant, dielectric properties and bulk modulus of bulk ceria calculations using IPs 1-6 and GGA+U are compared with experimental data in Table II.The lattice constant for IP1 and GGA+U are notably overestimated by > 1%, while for IP4 the lattice constant is underestimated by ∼0.5%.The remaining interatomic potentials reproduce the experimental lattice constant.Dielectric properties are reasonably reproduced by IPs 2-5, but are over-/under-estimated by IP1/IP6, respectively.All IPs overestimate the bulk modulus and are considerably larger compared to experimental data.Different lattice constants and bulk moduli are obtained for IP1 and IP2 due to the inclusion of the C ij term in IP2 as described above and the different parameters (Table I).As expected, the inclusion of the k 2 term in the shell model, illustrated by comparing IP2 and IP3, does not modify the structural, permittivity or bulk moduli properties of bulk ceria as the core and shell are superimposed in the bulk.However, altering the Buckingham potential parameters A ij , ρ ij (Eq.( 1)) and k 1 (Eq.( 2)) for only the Ce 4+ -O 2− interaction substantially influences the permittivity properties and bulk moduli, as shown by IP5 and IP6.All theoretically determined bulk properties calculated in this study are identical to those previously calculated and presented in the literature for each individual IP [10][11][12]19].

Surface Energies
The energies of the relaxed surfaces for IPs 1-6 and GGA+U are presented in Table II.For all methods, the relative stability of the surfaces, measured by the surface energy, decreases in the order (111) > (110) > (100).Physically this is what is expected, as the type III (100) surface undergoes significant reconstruction which will result in a greater energy loss upon surface formation and hence a higher surface energy.There is a considerable variation between the different IPs in the energies of the respective ceria surfaces.Differing only in the k 2 term of the shell parameter, IP2 and IP3 present identical surface energies.Comparing IP4 with IP5, despite displaying very different potential parameters, these IPs demonstrate almost identical surface energies for (111) and (110).IP5 and IP6 illustrate quite different surface energies, despite differing in only the A ij and ρ ij terms for the Ce 4+ -O 2− interaction and the k 1 shell model parameter.Comparing the different potentials, with the exception of IP1, IP6 gives considerably larger surface energies and IPs 2-3 illustrate consistently lower surface energies.The GGA+U surface energies for all ceria surfaces studied are consistently smaller than those calculated using IPs 2-6.However, IP1 gives rise to much lower surface energies compared to the other potentials and demonstrates some similarities to the GGA+U data (Table II).The (111) and   (110) surface energies of ceria have been calculated at the Hartree-Fock (HF) level by Gennard et al. [23], with a correlation correction of the Perdew-Wang functional added a posteriori.With the correlation correction the HF surface energies of ceria are very similar to IP6 described in this work.The correlation correction described by Gennard and coworkers increases the surface energy, however, they argue that the proper inclusion of correlation would be expected to lower the surface energies compared to IP data.This conclusion is upheld by Skorodumova et al. during the calculation of ceria surface energies using LDA+U /GGA+U where they also note smaller surface energy values compared to IP results [25].For other oxide materials, e.g.Al 2 O 3 , ZnO, TiO 2 and SnO 2 , DFT surface energies are also found to be smaller then HF and IPs [33].

Structural Relaxations
For the analysis of the surface relaxations, we highlight quantitative differences between the structural relaxations obtained with the IPs and GGA+U by considering the relaxation of the Ce and O ions perpendicular to the surface (Table III) and the contraction/elongation of the Ce-O bond distances (Table IV) in the surface and sub-surface layers of low index surfaces.This allows us to understand how well the structural relaxations computed with the IPs compare to those determined using GGA+U .Throughout the following discussion on the structural relaxation of ceria surfaces, the surface and sub-surface Ce and O ions are referred to as O1, Ce2, O3 and Ce4.This nomenclature is described in pictorial form in Fig. 1.
Table III illustrates the difference in core ion displacement perpendicular to the surface upon relaxation of the ceria surfaces with IPs 1-6, compared to our GGA+U results (similar results are observed for the shells).Negative values indicate a displacement into the bulk and positive values indicate that the ions relax out of the surface.As a whole, the magnitude of the ionic displace- ments perpendicular to the surface increases in the order (111) < (110) < (100), mirroring the respective surface energy and stability.Almost all IPs illustrate a displacement of the O ions that is larger than for the Ce ions, consistent with earlier IP [10-12, 15, 16], HF [24] studies and our GGA+U [28,29] results.For all three ceria surfaces, the ions of the surface and immediate sub-surface layers polarize, resulting in a separation of core and shell, which decays towards the centre of the slab.This displacement represents the polarization of the atom with the shell (valence electrons) displaced towards areas of positive charge and the core to areas of negative charge.The (100) surface undergoes the largest rearrangement of all the surfaces studied due to the removal of 50% of surface oxygen atoms to eliminate the dipole moment.Therefore, the displacement of the oxygen ions in ( 100  is notably larger than for the cerium ions.Owing to the availability of additional space on the surface due to the considerable surface re-arrangement in the (100), the second oxygen layer from the surface is split; one half of the ions shift outwards to the surface and the other half relaxes into the bulk.IPs 4-6 demonstrate that the ion displacement determined for all surfaces from the core positions are in good agreement with the GGA+U results.IPs 2-3 show almost identical displacement to that calculated via GGA+U in surface (111).However, these potentials overestimate Ce and underestimate O ion displacement in the other surfaces compared to the GGA+U results.IP1 causes some over-/under-estimation in ion displacement throughout all surfaces, but this is more pronounced for   IV), compared to other surfaces.The second Ce-O distance (Ce2-O3, Fig. 1, Table IV) remains unchanged in the GGA+U calculation, but is slightly relaxed for IPs 1-4 and slightly contracted for IPs 5-6.This overall small change is expected considering that this surface is generally considered the most stable and undergoes the least reconstruction.However, it should be noted that IP1 strongly overestimates the contraction/elongation for the first/second Ce-O bonds, respectively.IPs 2-3 and 5-6 separately produce identical surface and sub-surface bond distances, in comparison with each other, mimicking the trend described previously by their respective surface energies.The best agreement in the respective Ce-O distances for surface (111) is found between GGA+U and IP2 and IP3.A contraction of all Ce-O distances is observed for IPs and GGA+U , on surface (110), compared to bulk distances.Contraction of the second Ce-O distance (Ce2-O3, Fig. 1, Table IV) for IPs 5-6 is considerably overestimated.In general, the most comparable Ce-O distances to GGA+U are demonstrated by IP2 and IP4 for the surface (110).
Comparing the Ce-O distances of surface (100) to bulk data, GGA+U and all IPs produce a contraction of the first Ce-O bond, an elongation of the second and a further contraction of the third (O1-Ce2, Ce2-O3o and Ce2-O3v, respectively, Fig. 1, Table IV).The geometrical differences between Ce2-O3o and Ce2-O3v are illustrated in Fig. 1 and are determined by whether the oxygen sits below a surface oxygen (O3o) or below a space in the lattice (O3v) of the (100) surface.The contraction/elongation of the Ce-O distances is much larger in this surface, compared to both (111) and (110), due to the stronger surface reconstruction that occurs in (100) in order to remove the dipole moment.For this surface almost all of the Ce-O distances described with IPs are overestimated compared to GGA+U .As previously observed for the (111) surface, IPs 2-3 generate identical Ce-O distances compared to each other.The IP with the most compatibility to GGA+U and the least amount of under-/over-estimation of the data, for the (100) surface of ceria, is IP4.
B. Studying Ceria Surfaces using embedded-cluster QM/MM calculations Through our research testing the different IPs available in the literature (section III.A) against comparable GGA+U data for the low index ceria surfaces of (111), (110) and (100); we believe that of all the IPs studied here, the most comparable to our GGA+U data is IP4 (the potential of Grimes and Catlow [17]).While certain inconsistencies were established between the surface data modeled using these different computational methods, overall, IP4 demonstrated the most general consistency with GGA+U results.Based on this conclusion, we carried out trial QM/MM calculations on the (110) surface of ceria.Using IP4 to treat the MM section of the calculation, we investigated the stoichiometric (110) surface using a Ce 7 O 30 cluster and the (110) oxygen vacancy defect surface using a Ce 7 O 29 cluster.In this manner we were able to calculate the resulting defect energy for comparison with published periodic GGA+U [27,29] data.While Nolan et al. [29] published a p(2 × 1) defect calculation on the (110) surface of ceria with a corresponding defect energy of 1.99 eV, Fabris et al. [27] demonstrated a more isolated p(2 × 2) oxygen vacancy defect on the (110) surface to produce a defect energy of 1.57 eV.Our trial QM/MM calculations on the (110) surface using the Hartree-Fock (HF) approach yielded a completely isolated oxygen vacancy defect, as is the nature of the cluster calculation, with an accompanying defect energy of 1.05 eV.Overall, when considering the increasing isolation of the calculated defect, comparing the work of Nolan and Fabris with our own research, our QM/MM HF defect energy is comparable to results determined with GGA+U .
When we analyse the atomic spin population resulting from the defect QM calculation, we can readily observe two extra electrons, resulting from the surface removal of a neutral oxygen, on two adjacent cerium ions (Fig. 2)one on the surface and the other in the sub-surface layers of (110).Interestingly, these are not the same two reduced ions (Ce 3+ ) observed using GGA+U [29], and this is the first time modeling has identified subsurface reduction.This demonstrates that within the embedded cluster QM/MM calculation, the removal of an oxygen atom also models the simultaneous reduction of the surface (2Ce 4+ → 2Ce 3+ ).It should be noted that Herschend et al. [30] have carried out a similar HF embedded cluster calculation on (110) using a much simpler hemi-spherical design.Unlike our QM/MM model (described in section II.B), Herschend's sphere consists of a central cluster (Ce 10 O 20 ) of QM ions followed by a Ce interface region surrounded by point charges.This embedded cluster calculation produced a defect energy of 1.91 eV, indicating the importance of the long range polarization in defect energies.

IV. CONCLUSIONS
While the study of large systems can be performed with interatomic potentials (IPs), it is important to first validate their applicability to situations far removed from the bulk structures to which IPs are traditionally fitted.In this work we have examined how a number of IPs compare to an ab initio method for the description of the energetic and structural properties of bulk, and (111), ( 110) and (100) surfaces of ceria.Since GGA+U gives a good description of the structure and energetics of oxide surfaces in the literature, it was chosen as the ab initio method used for comparison.
Although both approaches give qualitatively similar results for surface stabilities and qualitative surface relaxations, there are quantitative differences.The GGA+U surface energies are smaller than IP and HF results, consistent with studies of other oxides.The structural properties of the IPs (1-6) studied are rather different and this is particularly true for the three low index surfaces of ceria studied.For all surfaces considered, as well as for bulk properties, the Butler derived potential (IP1) gives poorest agreement with the GGA+U results.The poor agreement of this potential with GGA+U and other IPs may be a consequence of the very large, positive shell charge in this potential.
We conclude that no single IP is able to describe the properties of pure ceria to a consistent degree of accuracy relative to the GGA+U benchmark.This can be traced back to the derivation of the IP, usually from fitting to bulk properties and may not be suitable in the asymmetric electric field at a surface.With this in mind, further investigation of the fitting of ceria interatomic potentials to accurate ab initio data is warranted.However, the best agreement is currently found between GGA+U and the interatomic potentials 2 and 4. While the surface energies of low index ceria are better represented by IP2, structure relaxations of the surfaces of ceria are modeled consistently more accurately, in comparison with GGA+U , when using IP4.
The periodic GGA+U method can accurately calculate http://www.sssj.org/ejssnt(J-Stage: http://www.jstage.jst.go.jp/browse/ejssnt/) e-Journal of Surface Science and Nanotechnology Volume 7 (2009) ceria defects on a periodic surface, however, due to the nature of the calculation these defects are usually on a crowded surface, i.e. a surface that contains many defects per surface area.Though in general IPs adequately model ceria surfaces and defect calculations, more accurate approaches are required in the future when considering isolated surface defects.In light of this we consider the more accurate embedded-cluster QM/MM calculations and have attempted trial calculations using a Ce 7 O 30 cluster to model the (110) surface of ceria.Based on these trial calculations we can compare the resulting defect energy with comparable published GGA+U results.Using this method we have established that the embedded-cluster QM/MM approach produces defect energies in line with GGA+U data and that this method also correctly models the subsequent reduction of the surface.
Indeed, our trial calculations indicate that subsurface reduction may be more energetically favourable than surface reduction.Future studies will investigate QM/MM embedded-cluster calculations on different surfaces of ceria utilizing several QM techniques and cluster sizes.

a
Short-range interaction cut-off distance 20.0 Å
the percentage difference in the Ce-O distances upon surface relaxation, compared to the respective bulk Ce-O distance.Negative values indicate a bond contraction and positive values indicate a bond relaxation.The (111) surface of ceria, compared to bulk distances, produces a small contraction of the first Ce-O surface bond (O1-Ce2, Fig. 1, Table

TABLE III :
Comparison of the difference in the core ion displacement ( Å) perpendicular to the surface, compared to GGA+U , of the interatomic potentials (IPs) 1-6 for Ce and O atoms in surface and first sub-surface layers.Positive difference in displacement indicates a relaxation of ions away from slab bulk-like centre and negative difference indicates a displacement of ions into the slab.

TABLE IV :
Comparison in the percentage (%) change in corecore Ce-O bond distances in the surface and first sub-surface layers of the (111), (110) and (100) surfaces compared to the bulk, calculated with interatomic potentials (IPs) 1-6 and GGA+U .A negative/positive % change indicates a contraction/elongation of the bond, respectively.