Self-Consistent-Charge Density-Functional Tight-Binding (SCC-DFTB

Subscriber access provided by University of Newcastle, Australia

Article

Self-Consistent-Charge Density-Functional TightBinding (SCC-DFTB) Parameters for Ceria in 0D to 3D Jolla Per Kullgren, Matthew Jason Wolf, Kersti Hermansson, Christof Köhler, Bálint Aradi, Thomas Frauenheim, and Peter Broqvist J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.6b10557 • Publication Date (Web): 02 Jan 2017 Downloaded from http://pubs.acs.org on January 3, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry C is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Self-Consistent-Charge Density-Functional Tight-Binding (SCC-DFTB) Parameters for Ceria in 0D to 3D Jolla Kullgren,∗,† Matthew J. Wolf,† Kersti Hermansson,† Christof Köhler,‡ Bálint Aradi,‡ Thomas Frauenheim,‡ and Peter Broqvist∗,† Department of Chemistry–˚ Angström, Uppsala University, Box 538, S-751 21, Uppsala, Sweden., and Bremen Center for Computational Materials Science, Universit¨ at Bremen, P.o.B. 330440, D-28334 Bremen, Germany. E-mail: [email protected]; [email protected]

∗

To whom correspondence should be addressed Department of Chemistry–˚ Angstr¨ om, Uppsala University, Box 538, S-751 21, Uppsala, Sweden. ‡ Bremen Center for Computational Materials Science, Universität Bremen, P.o.B. 330440, D-28334 Bremen, Germany. †

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract Reducible oxides such as CeO2 are challanging to describe with standard density functional theory (DFT) due to the mixed valence states of the cations, and often require the use of additional correction schemes, and/or more computationally expensive methods. This adds a new layer of complexity when it comes to the generation of Slater-Koster tables and the corresponding repulsive potentials for self-consistent density functional based tight binding (SCC-DFTB) calculations of such materials. In this work, we provide guidelines for how to set up a parameterisation scheme for mixed valence oxides within the SCC-DFTB framework, with a focus on reproducing structural and electronic properties as well as redox reaction energies calculated using a reference DFT method. This parameterisation procedure has been used to generate parameters for Ce–O interactions, with Ce in its +III or +IV formal oxidation states. The generated parameter set is validated through comparison to DFT calculations for various ceria (CeO2 ) and reduced ceria (CeO2−x ) systems of different dimensionalities ranging from 0D (nano-particles) to 3D (bulk). As oxygen vacancy defects in ceria are of crucial importance to many technological applications, special focus is directed towards the capability of describing such defects accurately.

Introduction Ceria (CeO2 ), is used in many different technological applications, but is probably most well known for its important role in the three-way catalyst used to clean exhaust fumes from automobile engines 1 and as a solid-state electrolyte in solid oxide fuel cells. 2 In these applications, ceria is used as an oxygen buffer material, where the most relevant processes mainly rely on ceria’s ability to release and store oxygen by the creation and subsequent healing of oxygen vacancies within the structure. The oxygen buffering is made possible by the ease with which Ce can toggle between the +IV and +III formal oxidation states i.e between Ce4+ and Ce3+ . Consequently, oxygen vacancies in ceria have been intensively

2


Page 2 of 55

Page 3 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


studied, both experimentally and theoretically, in the past few decades (See for example the review by Paier et al. 3 ). Recent advances in nanotechnology allow nowadays for the synthesis of ceria nanoparticles of controlled sizes and shapes. 4–9 Studies of, for example, the oxygen storage capacity for such nanoparticles have shown an enhancement of this property for particles with diameters less than 5 nm, which has been explained to be due to their ability to transform adsorbed 10–12 oxygen molecules to superoxide ions (O− Those studies suggested that for small sizes, 2 ).

ceria becomes redox active at lower temperatures and without the need to form explicit oxygen vacancies. New nanoscale findings may expand the use of ceria beyond “standard” ceria applications. One such example is the potential use in nanomedical applications where ceria nanoparticles have been demonstrated to act as therapeutic agents in the treatment of oxidative stress in living cells. 13,14 Despite the intense interest in ceria based materials, theoretical methods capable of modelling the chemistry of realistic systems, both in terms of size and time scales, are scarse. Although classical force field methods (FFs) (e.g the potentials presented in Refs. 15–21 ) can be used to simulate relatively large and structurally complex ceria nanoparticles (see e.g. Ref. ? ), they are not able to describe the breaking and formation of chemical bonds, nor changes in the oxidation state of ions, both of which are central effects in redox chemistry. Here, more advanced FF formulations, which allow for charge transfer and bond breaking (socalled reactive force fields), open up new possibilities. One such formulation is ReaxFF, 22–24 which was recently used to model ceria and reduced ceria. 25 However, describing localised electrons associated with Ce3+ ions using this family of FFs is challenging, and consequently also the detailed structure around vacancies is not fully captured. In general, quantum chemical methods such as Density Functional Theory (DFT), are of course able to describe such effects very well, but due to their relatively high computational cost, they are only feasible to use for quite small, and not too complex, ceria systems. Thus, there is a need for approximate methods to enable studies of large systems which still allow for an accurate

3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

description of charge transfer reactions. One such method is the so-called Self Consistent Charge Density Functional Tight Binding method (SCC-DFTB) which is an approximate, parametrised form of DFT. 26,27 The method is based on a second-order expansion of the total energy in terms of deviations from a reference electronic density, which is typically taken to correspond to a superposition of atomic densities. Hamiltonian matrix elements describing the electronic structure of the system are pre-calculated using DFT for each distinct pair of atomic species in the system, and stored in so-called Slater-Koster (SK) tables. The total energy and forces are then corrected with a two-body force-field term, the so-called repulsive potential, which is typically fitted so as to give the best compromise with respect to the description of key properties (typically energies and/or structures) of a set of systems calculated using DFT. Owing to the fact that the matrix elements are not calculated on the fly, SCC-DFTB is at least two orders of magnitude faster than standard semi-local DFT. Published SCC-DFTB parameter sets for oxide materials often yeild cohesive energies (atomisation energies) that are severely overestimated, 28,29 while the binding energies of simple molecules (such as the O2 molecule) are well described when compared to reference DFT data. 26 Such inbalance will lead to errors compared to DFT in e.g vacancy formation energies or reduction energies, which often are central quantities in redox chemistry. At the same time we note that such properties are allready difficult to model using (semi-)local DFT, and often requires extensions of DFT, such as the DFT+U method, or higher levels of theory altogether. Thus, capturing such properties within the SCC-DFTB formalism adds a new layer of complexity when it comes to the parameterisation, and necessitates the development of new approaches for parameterisation to enable computationally efficient and accurate studies of redox behaviours in this class of materials. In this work, we present parameters suited for studying redox properties of ceria using SCC-DFTB, and discuss the approach that we used to generate them. The performance of the parameters is validated through comparisons to DFT by systematically testing using

4


Page 4 of 55

Page 5 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


a number of CeO2 structures of different complexities and stoichiometries (bulk, surfaces, surface steps and nanoparticles). It is our hope that the parameterisation effort provided here in will provide guidelines for how to generate SCC-DFTB parameters for mixed valence, reducible oxides in general, and constitute a first step towards an automatic parametrisation procedure for such materials. This paper is organised as follows. In Sec. , we present the computational details and theoretical methods used. Sec. discusses the parameterisation approach we have undertaken to generate the SCC-DFTB parameters. In Sec. , we thoroughly test the generated parameters for stoichiometric and partially reduced ceria systems of different dimensionality (0D to 3D). Sec. provides some concluding remarks.

Computational details In this section, we will describe the theoretical methods, both the reference DFT method and the SCC-DFTB method, that we have generated.

Density Functional Theory The reference method, to which we compare all our results, is density functional theory in the implementation with plane waves and pseudopotentials. We used the GGA–PBE functional developed by Perdew et al. 30 The pseudopotentials were generated in the form of the Projector Augmented Wave (PAW) method described by Blöchl 31 to treat the core electrons. The valence space consists of 6 and 12 electrons for oxygen and cerium, respectively. In addition, to accurately describe the strongly correlated Ce4f states, we have used the DFT+U method from Dudarev et al. 32 Motivated by previous work, a U -value of 5 eV was used which gives a good overall description of ceria. 33–37 All calculations were done using the VASP software package. 38–41

5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 55

SCC-DFTB The SCC-DFTB total energy is derived from a second-order expansion of the DFT total energy expression with respect to charge density fluctuations. 27 The total energy can be expressed as a sum of three terms according to:

ESCC−DF T B = EBS + E2nd + Erep .

(1)

EBS is the band structure term and is obtained from the sum of all occupied eigenstates. The second order term, E2nd , is the energy caused by density fluctuations (deviations from the reference atomic charge densities). The final term, Erep , is dominated by the ion–ion repulsion, hence the name repulsive energy, but it includes all remaining interactions within the DFT total energy. In practice, Erep is expressed as a sum of pair-wise 2-body terms which we here denote the repulsive potential Vrep (r). In the SCC-DFTB formalism, the electronic wavefunctions used in the calculation of EBS are expressed in a minimal basis of pseudo-atomic orbitals generated from an all-electron DFT calculation for each atom. The distance dependent overlap matrix elements and the matrix elements representing the Hamiltonian in this basis are pre-calculated once and stored in so-called Slater-Koster tables. The second-order term, E2nd , makes use of spherically symmetric, exponentially decaying atomic charge densities, the magnitudes of which are determined by the respective Mulliken charges, and which determine the on-site and inter-atomic contributions to the second-order term. These contributions are expressed in terms of an l-dependent internal Hubbard U parameter, being twice the atomic absolute hardness, which is also tabulated along with the Slater-Koster tables. The Slater-Koster tables and repulsive potentials Vrep (r) are unique to each atom–atom pair. For ceria, we have three unique pairs corresponding to Ce–Ce, Ce–O and O–O. Out of these three, only parameters for the O–O interaction have been presented in the literature.

6


Page 7 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Here, we have used the O–O parameters from the mio-1-1 set 26 and have developed new parameters for Ce–Ce and Ce–O. For strongly correlated electrons, such as Ce4f electrons, SCC-DFTB suffers from the same problems as its parent DFT method does, namely the tendency to over-delocalise such states. This implies that we also need a +U correction for SCC-DFTB calculations. In this work, we have used the fully localised limit form of the +U correction, as implemented by Hourahine et al. 42 All SCC-DFTB calculations were performed using the DFTB+ code. 43

SCC-DFTB parameterisation for ceria in this project When using an approximate, parameterised method such as SCC-DFTB, one would ideally like to have an automatic, general purpose procedure for generating the parameters that leads to an optimal set for any system. While much progress has been made towards achieving that goal, 44–46 there remain difficult cases for which one must take a more hands on approach, and to a greater extent rely on chemical intuition and experience. In general, a parameterisation for any atom pair is started by performing all-electron reference calculations of the atoms in their ground states. The atomic references are then straightforwardly used to generate SK-tables which reproduce the electronic properties of a reference material or molecule of choice calculated with the parent DFT method. In the case of metal oxides this reference material is typically the corresponding bulk metal and the bulk oxide. Thereafter, a repulsive potential Vrep (r) is fitted to ensure that the structural properties of some reference structures agree. This parameterisation procedure works well in most cases, and in the density superposition scheme 26,27 which is used in current publication it involves, in principle, two parameters for the Hamiltonian, the compression radii for the electronic wave functions (rΨ ) and the compression radii for the density (rn ). The compressed density is used to calculate the poten-

7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

tials in the overlap and Hamiltonian integrals. Additionally, one has the repulsive potential Vrep (r) which corrects the total energy. However, when using this standard procedure for the parameterisation of Ce–O interactions we immediately ran into problems, particularly when it came to reproducing reduction energies for ceria systems of different dimensionalities. Therefore, we found it necessary to modify the parameterisation procedure following a more pragmatic approach, in which we first set up a list of target properties, which we wanted to be able to reproduce with the desired level of accuracy, followed by a systematic search with an extended set of adjustable parameters. We will here further discuss the parameterisation approach we have followed in more detail starting with the targets, i.e. the properties we want to be able to model as accurately as possible. In the second section we will discuss and motivate our choice of all adjustable parameters. Finally, we will describe the process of obtaining the optimal parameters for simulations of the redox behaviour of ceria.

Target properties Before starting the parameterisation, we set up a list of ceria systems and target properties, similar to a training set commonly used in the parameterisation of force fields, that we want to be able to study using the SCC-DFTB method, namely the following basic and crucial proporties which are of course also essential for the technological applications discussed in the Introduction: • Relative stabilities of two bulk CeO2 polymorphs with different coordinations, namely the rutile and fluorite structures. • Surface energies of the low index surfaces of ceria in the fluorite structure. • Vacancy formation energies in bulk ceria in the fluorite structure and at its low-index surfaces. • Stability of fluorite ceria nano structures. 8


Page 8 of 55

Page 9 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


• Electronic structure, primarily the electronic Density of States (DOS) diagram, both of stoichiometric systems and of systems with oxygen vacancies, namely partially reduced ceria systems (here denoted CeO2−x and the fully reduced Ce2 O3 , Regarding the bulk structures, we note that, while the existence of a rutile structured bulk phase has not been observed experimentally, it has been speculated, based on theoretical calculations that nanorods can adopt this structure. In any case, we use the rutlie structure with Ce ions of lower 6-coordination (6C) vs. 8-coordination (8C) for the fluorite structure, since we found to be important in other parameterisations, e.g. of ZnO, 47 to include a range of coordination numbers. Note also that the Ce metal, which would normally be an important target in a typical SCC-DFTB parameterisation, has been omitted from the above set since we found at an early stage that it was problematic to find suitable parameters that achieved adequate transferability with respect to all of the oxidation states of Ce, from 0 to +IV. All the entries in the list above were treated as being equally important at the beginning of the parameterisation. However, during the course of parameterisation we needed to alter this, and at the final stages of parameterisation, we put extra emphasis on properties related to the fluorite phase of ceria (entries 2 to 5 in the list above).

Adjustable parameters In the standard cases of SCC-DFTB parameterisation, the adjustable parameters can be divided into two categories. First, there are the parameters that directly affect the tabulated integrals and other data stored in the Slater-Koster tables. These parameters are, normally, the compression radius for the wave functions (rΨ ) and the density (rn ). In our parameterisation, we have used quadratic compression potentials and relativistic effects were taken into account according to the so-called ZORA scheme. Secondly, we have parameters describing the repulsive potential. To allow for flexibility while still retaining a smooth repulsive potential, we describe the repulsive potential with 9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

an analytical function that was fitted using the General Utility Lattice Program (GULP). 48 We use the four range Buckingham expression as implemented in GULP and follow the procedure outlined in Ref., 47 where this form for the repulsive potential was shown to be powerful and work well for a qualitative description of different phases of ZnO. 47 Moreover, we have found it necessary to also allow for the modification of other parameters which concerns the generation of the Slater-Koster tables which are normally obtained from all-electron calculations for the respective atom in its ground electronic state. Here, we found that changing the reference electronic state of the Ce atom improved the description of our calculated reduction energies and we optimised the Ce electronic configuration by allowing charge transfer between the Ce5d and Ce4f states and by allowing for a fractional occupation of the electronic states of the Ce atom. In practice, we found that this modification mainly affects the tabulated atomic eigenstate energies and the atomic absolute hardness. To justify this procedure, we note that the description of the electronic configurations of the cerium atom is a challenging task even for post-Hartree-Fock methods. In post-HartreeFock theory the experimental d-f charge transfer energy, i.e. the energy difference between atomic f n d1 s1 and f (n+1) d0 s2 configurations, can be recovered, while in relativistic DFT calculations the energy difference has the wrong positive sign, thereby favouring the f 2 d0 s2 cerium configuration. 49,50 This situation is reminiscent of the s-d charge transfer energies in the 3d transition metal atom series. 51–53 Here Painter et al. 54 actually found that the atomic ground-state of Fe in LDA is non-integer and corresponds to a 3d6.384 s1.62 state. Similarly for the Ce atom using PBE we find a [Xe]4f1.6 5d0.4 6s2 ground-state, rather than the experimentally suggested [Xe]4f1 5d1 6s2 one. Thus, motivated by this finding we decided to allow for a non-integer occupation number when generating our Slater-Koster tables. Lastly, as discussed in the Method section, localised f -electrons associated with Ce3+ ions are reduced ceria is not correctly described niether using DFT within the local density (LDA) nor with the general gradient approximation (GGA). The failure is caused by a

10


Page 10 of 55

Page 11 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


large self-interaction error (SI-error) for the spatially confined f -electrons which leads to unphysical delocalisation. SCC-DFTB based on GGA will suffer from the same problems. One conceptually simple, and computationally efficient, way to correct for the SI-error is the so-called Hubbard (+U ) correction. In this approach, the f -electrons are projected out and treated in a Hubbard-like model with an on-site Coulomb interaction controlled by the parameter U . This method imposes an energy penalty on partial f -occupancy and drives the DFT ground state away from the delocalised LDA/GGA ground state, towards an electronic structure that more closely corresponds to an integer f -occupancy. In the current work we include the +U correction in both the SCC-DFTB simulations as well as in the reference DFT calculations. We note that while, qualitatively, the +U correction in SCC-DFTB has the same effect as in DFT, i.e. it localises the Ce(4f ) electron, shifts the associated eigenvalues and leads to an expansion of the lattice, the effect of the external +U parameter is found to be larger in the SCC-DFTB calculations (given a fixed repulsive potential), both on the shifted eigenvalues and the lattice expansion. Consequently, the +U value used in the SCC-DFTB was set to be lower than that used in the DFT calculations. Because of these differences, we were forced to generate a new repulsive potential for each U-value to reproduce the reference data. Thus, the +U correction was added as a parameter that has to be optimised along with the other considered parameters. We note here that this means that the +U -value cannot be used to tune vacancy formation energies a posteriori once all other target properties have been optimized. This is not least due to the fact that a large part of the vacancy formation energy originate from structural relaxation which is indirectly affected by the +U -value. In addition, we found it necessary to also use a +U on the Ce5d states to obtain electronic results in better accordance with the reference data. In this work, we did not tune the +U value for the d-states, but instead used the same value as for the Ce4f states. In summary, the adjustable parameters in the parameterisation of ceria are the following: • d-f charge transfer in the Ce atom (focc − docc ) 11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

• Compression radii for atomic wavefunctions (rΨ ) • Compression radii for the density (rn ) • External Hubbard U parameter for Ce4f and 5d (U ) • Repulsive potential between Ce-Ce and Ce-O (Vrep (r)) Finally, during the course of the parameterisation, we made the following observations regarding how the different parameters affect properties of the Ce-O system. We found that the compression radii had the largest affect on the band structure of bulk CeO2 while surface energies and nanoparticle stabilities were most strongly affected by the repulsive potential. For reduced ceria, the band structure and energetics was found to be dependent also on the external +U parameter, the atomic eigenstates and the atom absolute hardness. As such, the adjustable parameters can, to some extent, be viewed as belonging to quasi-independent sets and be tuned individually to optimise certain properties. However, for a simultaneously accurate description of both ceria and reduced ceria systems of various dimensionality, the parameterisation becomes highly non-linear and it becomes indeed challenging to obtain a balance between the different adjustable parameters. We believe that we have met this challange rather well.

Finding the optimal parameters In our inital attempts at parameterisation, we found it rather straightforward to generate sets which were able to reproduce, e.g., the band structure of ceria, or the bulk vacancy formation energy, but often at the expense of the description of other properties, such as the surface energy. To produce a set of parameters that provides a well balanced description of the broad range of ceria properties under consideration here, we used an iterative parameterisation procedure as illustrated in Fig. 1. In this procedure, we start by optimising the different compression radii to reproduce band structures of some test cases with fixed geometries, typically stoichiometric bulk ceria 12


Page 12 of 55

Page 13 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and bulk ceria with an oxygen vacancy. When a sufficient quality was reached, we continued to fit repulsive potentials, Vrep (r), for Ce–O and Ce–Ce, respectively. For the Ce–O repulsive potential, we used volume scans for the two phases of CeO2 targeted. As mentioned we did not include the Ce metal in our fitting, but instead we used a bond scan of the Ce2 O3 molecule to fit the repulsive potential for Ce-Ce. After this step, we continued into the testing phase where all the other ceria systems and properties in our target list were scrutinized. The whole procedure was then repeated until satisfactory results for all desired properties was obtained. We found that it is important to include full relaxation in energy comparisons between SCC-DFTB and DFT of, e.g. vacancy formation and surface energies. Such relaxation cannot be done before having a repulsive potential in place. This prohibited us at this point from using an automatic fitting procedure in this work. The aggressive tuning of the parameters to reproduce outlined target properties with acceptable accuracy unfortunately leads to losses in transferability. We have, for example, intentionally sacrificed the description of the Ce metal during this process. Also, we had to sacrifice the target regarding the two polymorphs of ceria, since this was in conflict with the possibility of obtaining accurate surface energies. This is further discussed in the validation section below. Following the outlined iterative strategy (c.f. Fig. 1), we reached the following valCe(s,p,d)

ues for the parameters listed above: rΨ

Ce(f )

= 4.5 bohr, rΨ

=4.0 bohr, rVCe =12.0 bohr,

focc =1.35e− , docc =0.65e− and U =2.5 eV for both f and d states. Since the on-site Hubbard correction is crucial in our scheme, we hereon use the short term DFTB+U when referring to SCC-DFTB+U calculations.

13



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Validation of parameters Here, we demonstrate the applicability of the generated DFTB+U parameters for ceria by comparison of our SCC-DFTB results with the corresponding result from our reference method, PBE+U , for a number of important ceria systems. Since much of the chemistry of ceria crucially depends on the presence of oxygen vacancies as mentioned, we placed extra focus on the DFTB+U description of these defects during the parametrisation and on the transferability between models of different dimensionality (3D bulk to 0D nanoparticles). Here we will highlight these results.

Bulk phases – 3D In the validation of bulk properties, we considered both stoichiometric CeO2 , partially reduced CeO2−x , and fully reduced Ce2 O3 . Stoichiometric CeO2 In the cubic fluorite structure (Fig. 2a), which is the most stable structure for CeO2 under normal conditions, the oxygen ions are four coordinated and cerium ions are eight coordinated. The rutile structured phase (Fig. 2b), on the other hand serves here as a second reference because it exhibits coordination numbers of three and six for oxygen and cerium ions, respectively. Fig. 3 shows a DFTB+U /PBE+U comparison for energy-volume scans of the fluorite and rutile phases of CeO2 . Our parameters are seen to correctly predict the fluorite phase to be the most stable polymorph over a wide range of volumes. However, the quantitative agreement between the DFTB+U and PBE+U calculations is off with DFTB+U predicting a too stable rutile phase. While quantitative agreement would be desirable, we found that achieving this always worsens the description of the fluorite surfaces and nanoparticles. Thus, there seems to be a ”built in” favorisation of the rutile stucture in our DFTB+U approach.

14


Page 14 of 55

Page 15 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


We have not been able to fully resolve and understand the origin of this effect. In the Suporting Information to this paper we give a full comparisons of the current parameter set to a parameter set where the rutile structure is better described. There we discuss the trade-off in more detail and allow the reader to form his or her opionon of the problem. From a more practical point, the rutile structure is of lesser importance for the intended applications and we judge that it is necessary and sufficient to simply ascertain that this structure is significantly destabilised compared to the fluorite phase. We note that the fluorite lattice parameter obtained in the DFTB+U calculation is in very good agreement with the PBE+U data, as should be, given that it was used in the training of the parameters. For the same reason, the bulk modulues, calculated by fitting the data of Fig. 3 to the Birch-Murnaghan equation of state is also in very good agreement. Both of these data are given in Table 1. A schematic illustration of the electronic structure of ceria is shown in Fig. 4(a). The energy gap values indicated in the figure were taken from Ref., 34 where the available experimental measurements were scrutinised and summarised. The energy gap values ∆ǫ from the valence band maximum (VBM, being mainly composed of O2p states) to the localised Ce4f and Ce5d states, calculated using DFTB+U and PBE+U are given in Table 1. Considered togheter, the various gaps calculated using the DFTB+U method are in quite good agreement with the PBE+U data. Again, this is partly by construction, as these PBE+U quantities were taken into account when the SCC-DFTB parameters were generated. Compared to the experiments presented in Fig. 4(a), we note that both methods underestimate the experimental values. The oxygen vacancy in bulk ceria Oxygen vacancy formation in CeO2 , forming CeO2−x , involves the reduction of two cerium ions from Ce4+ to Ce3+ . This reduction reaction can be written as:

15



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 55

1 CeO2 (s) → CeO2−x (s) + O2 (g) 2

(2)

1 3+ 2− 2Ce4+ latt + Olatt → 2Celatt + O2 (g). 2

(3)

or, equivalently,

At Ce3+ sites, the extra electron occupies an otherwise empty Ce4f -state and the Ce4f band in ceria splits into occupied and unoccupied bands, which we label Ce4f (occ) and Ce4f (unocc) band. The two bands lie about 2 eV and 3 eV above the valence band, respectively, as shown schematically in Fig. 4(b). We have modelled a bulk oxygen vacancy, V, by removing one oxygen ion from a 2×2×2 crystallographic supercell (96 atoms) of ceria. The different structures were then optimised using both the DFTB+U and PBE+U methods. As has been the subject of much discussion in the literature, the formation energy of an oxygen vacancy depends on whether the Ce3+ ions are in nearest neighbour (NN) or next nearest neighbour (NNN) positions with respect to the vacancy (see Fig. 5(a) and 5(b)). The calculated vacancy formation energies are 2.93 (3.26) and 3.04 (3.22) eV using DFTB+U (PBE+U ) for VN N and VN N N , respectively, and are given in Table 2. The absolute values for the vacancy formation energies are of the right magnitude, but the errors with respect to our reference method are in the order of 0.2 to 0.3 eV and also the stability ordering between VN N and VN N N is at variance. The energy difference for the different localisations is, however, small (∼0.1 eV), for both methods. The formation of an oxygen vacancy and its accompanying Ce3+ ions is associated with a substantial relaxation of the local structure. We have quantified this relaxation by measuring interatomic distances between the four Ce ions closest to the vacancy. Figs. 5a and 5b show these four Ce ions and their relation to the vacancy. The distances are given in Table 3; the DFTB+U values are in very good agreement with those from PBE+U , with errors smaller than 1% on average.

16


Page 17 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Fig. 5 shows the electronic density of states (DOS) for VN N and VN N N calculated with the DFTB+U and PBE+U methods. The band widths obtained with the DFTB+U method are smaller than for the PBE+U method. A similar observation was made for ZnO in Ref., 47 and is due to the use of a minimal basis set in the DFTB+U method. Otherwise, as noted above, the electronic gaps between the valence band (consisting of mainly O2p states) and the localised empty Ce4f states are well reproduced with the DFTB+U method. However, while the filled Ce4f states fall at sligthly different positions in the band gap for VN N and VN N N with the PBE+U method, they are similar, within 0.01 eV with the DFTB+U method. The lower value obtained for VN N N is not reproduced with DFTB+U . Thus, comparing ∆ǫ for O2p to occupied Ce4f , DFTB+U slightly overestimates the PBE+U calculated difference by 0.47 eV (0.25 eV) for VN N (VN N N ), respectively. Bulk Ce2 O3 Most technical applications of ceria make use of its redox properties i.e its ability to toggle reversibly between Ce4+ and Ce3+ . The end product in a full reduction cycle would be Ce2 O3 , according to the reaction 1 2CeO2 (s) → Ce2 O3 (s) + O2 (g). 2

(4)

We have therefore tested our new DFTB+U parameters also for this reaction. Two known phases of Ce2 O3 have been discussed in the literature, the A- and the C-types, where the A-Ce2 O3 is typically observed in experimental studies and is the most stable. We have only considered the A-Ce2 O3 in this study, the structure of which is shown in Figure 2c. The calculated lattice constants a and c are 3.82 ˚ A (3.92 ˚ A) and 5.98 ˚ A (6.18 ˚ A) using DFTB+U (PBE+U ). Both are in fair agreement with the experimental values of a=3.89 ˚ A and c=6.06 ˚ A. 55 The calculated energy of formation for reaction 4 is given in Table 2. There is a discrepancy between the DFTB+U and the reference PBE+U method. But both

17



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

methods, however, underestimate the experimental value of ∼4 eV. 56 The calculated electron densities of states using the two methods are shown in Figure 6. The agreement between DFTB+U and PBE+U is as good as that we saw in the previous DOS plots. The unoccupied d-states in DFTB+U are however shifted to higher values compared to PBE+U (not shown).

Ceria surfaces – 2D This section concerns the stabilities and structures of the low-index surfaces of ceria in the fluorite structured crystalline phase. As in the case of the bulk, we also consider partial reduction of the stoichiometric CeO2 systems by studying oxygen vacancy formation on these surfaces. Stoichiometric surfaces Three different low-index surfaces have been tested: (111), (110) and a reconstructed (100) surface (see Fig. 7). The bulk terminated (100) surface is polar and various compensation mechanisms to quench this polarity have been proposed in the literature. This will not be addressed here; instead we focus on whether or not DFTB+U can reproduce the results from a PBE+U reference structure. Therefore, we adopt a simple model for this surface consisting of ordered oxygen vacancies, and constructed from the bulk terminated surface by moving every other O ion row to the opposite face; it is henceforth referred to as (100)r. The surfaces were modelled using slabs of 7 triple atomic layers for the (111) and (110) surfaces, and by 11 atomic layers for the (100)r surface; see Fig. 7. In all models, a vacuum gap of at least 12 ˚ A was used to eliminate artificial interactions between periodic images along the z-direction. The structures of all surface models were optimised using both the DFTB+U and PBE+U methods, with the surfaces built from and kept at their respective optimised bulk lattice parameters. First, we compare surface energies Esurf calculated from: 18


Page 18 of 55

Page 19 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Esurf =

Eslab − N · Ebulk , 2A

(5)

where Eslab and Ebulk are the total energies of the slab and the bulk unit cell, respectively. A is the area of the exposed surface and the factor of 2 accounts for the fact that each slab exposes two symmetrically equivalent surfaces, one on the top and one on the bottom. The calculated surface energies for all three surfaces are given in Table 4. The agreement between DFTB+U and PBE+U surface energy for the (111) surface is overestimated with our DFTB+U parameters, but is overall very good. This overestimation is, as discussed above, is at least partly, a result of our efforts to ensure that the bulk fluorite structure is more stable than the rutile; aiming for a lower surface energy for the (111) surface was found to simultaneously result in a too stable rutile structure. The reported value in Table 4 is thus a result of a compromise where the description of the surface was emphasized. Importantly, the DFTB+U method is seen to well reproduce the PBE+U trends in the surface energies, with the (111) surface being most stable. The surfaces strucutures have been quantified by calculating the relaxation ρ and the rumpling ǫ. These are commonly defined as:

ρ=

1 danion + dcation − 2dbulk 2 dbulk

(6)

danion − dcation . dbulk

(7)

ǫ=

Note that here, for the (111) surface, danion is the distance between the O ions in the first atomic triple-layer and the corresponding ones in the second triple-layer (cf. Fig. 7), i.e. the difference in the z-coordinates of the O ions in atomic layers 1 and 4. Correspondingly, dcation is the difference in the z-coordinates of the Ce ions in atomic layers 2 and 5. The quantity dbulk is the corresponding geometry-optimised bulk Ce-O inter-layer distance in the direction normal to the surface. The calculated values are given in . The structural relaxation and 19



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

rumpling values predicted by the two methods are in very good agreement (Table 4). The two more open structures, (110) and (100)r shows slightly larger errors that (111). In summary, we find that with our generated DFTB parameters, we are able to overall reproduce the PBE+U surface energy well, the largest discrepency (0.93 vs 0.71 eV) being that of the most stable (111) surface. The structural relaxation and rumpling of all surfaces are very well reproduced. Surface and sub-surface oxygen vacancies at the CeO2 (111) surface The nature of surface vacancies (SV) and sub-subsurface vacancies (SSV) at the (111) surface has been studied extensively in the literature, using both experimental techniques and theoretical simulations. There are, however, still some unresolved discrepancies between the results of the experimental and theoretical studies, especially related to the relative stability of SV and SSV and their propensity to form clusters at the surface. In the current investigation we will not attempt to resolve these controversies but rather focus on a comparison between the DFTB+U and PBE+U descriptions of these defects by calculating their formation energies at the CeO2 (111) surface with different electron localisation patterns. As in the case of the bulk oxygen vacancy, previous DFT studies suggest that the vacancy formation energy depends on the locations of the Ce3+ ions, and that the most stable configurations are those with Ce3+ ions occupying next-nearest neighbour positions with respect to the vacancy. In this study we consider two configurations for each vacancy type, either with both Ce3+ ions in nearest (NN) or next-nearest neighbour (NNN) positions. In both configurations, shown in Fig. 8, the Ce3+ are located in the uppermost triple layer of the slab. The calculated vacancy formation energies are reported in Table 2. Clearly, DFTB+U is able to reproduce the vacancy formation energies of the PBE+U reference well, which has been one of the main goals for the parameterisation. Thus, both methods predict that SSV with NNN localisation is the most stable type of vacancy at the (111) surface. However,

20


Page 20 of 55

Page 21 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


for the other structures explored, DFTB+U predicts a similar formation energy of ∼2.4 eV, while PBE+U shows some variation (∼2.3-2.6 eV). We further analyse the local structure in the neighbourhood of the SV and SSV defects at the (111) surface. The PBE+U method show that the relaxation around the SV is characterised by a strongly asymmetric distortion of the surface oxygen ions surrounding the vacancy, leading to unequal distances between these pairs of ions. For the SSV, the relaxation is characterised by a protrusion of surface oxygen ions at next nearest neighbour position with respect to the vacancy. We have measured such pair-wise distances from the optimised geometries using both the DFTB+U and PBE+U methods as indicated in Fig. 8(a) and 8(b) for the SV and SSV, respectively. The results for the geometrical parameters indicated in the figure are given in Table 5 which shows that there is very good agreement between the two methods, as was also found for the bulk vacancy. For the SSV defect, the nearest neighbours are always at “higher altitude” than the next-nearest neighbours, but the difference is larger for the NNN configuration. This holds for both DFTB+U and PBE+U . We conclude that DFTB+U gives a good description of both the relative stability of and local structure around vacancies at the CeO2 (111) surface, i.e. it reproduces the PBE+U results well. For all optimised structures, we have additionally examined the electron DOS and the projected DOS on the Ce3+ ions. In all cases, the agreement between the two methods is good. We illustrate this for the most stable vacancy defect, namely the SSVN N N , in Fig. 8(c). The surface oxygen vacancy at the CeO2 (110) surface The oxygen vacancy at the (110) surface has attracted considerably less attention than the one on the (111) surface. This is likely related to the difficulties in obtaining well ordered (110) surfaces experimentally and to the fact that the (110) surface normally represents a minority facet. However, a number of DFT investigations have carefully characterised the

21



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

atomic and electronic structure of surface oxygen vacancies at the (110) surface. 57–59 Two different types of vacancies have been identified. One is essentially just the pristine surface with one oxygen ion missing (in-plane vacancy, SVIP V ) while the other structure involves a large displacement of a neighbouring oxygen ion which occupies a bridging position (bridge type vacancy, SVBT V ). There is a strong interplay between the atomic relaxation pattern and the locations of the two excess electrons accompanying the vacancy, often called the electron localisation pattern. In the current investigation we compare two such patterns, namely the one with lowest energy in Ref. 59 which is of the SVBT V type and the structure in which both Ce3+ are nearest neighbours to the vacancy, which is of SVIP V type. These two structures are shown schematically in Fig. 9(a) and 9(b), respectively. The calculated formation energies are reported in Table 2 for both methods. They both predict the formation energy on the (110) surface to be lower than that for the vacancies on the (111) surface. Overall, the agreement between DFTB+U and PBE+U is (only) fair for the (110) surface. We have further analysed the local relaxation around the SVIP V and SVBT V at the (110) surface in terms of the displacement of the closest oxygen ion (see Fig. 9). We have measured the distance when projected onto the (110) surface plane. For the SVIP V and SVBT V structures, we obtain values of 0.22 (0.24) and 1.36 (1.37) ˚ A, respectively, using the DFTB+U (PBE+U ) method. Thus, the local relaxation induced by the vacancy formation at the (110) surface is well reproduced by our DFTB+U parameters. We have also analysed the electronic total DOS and DOS projected on the Ce4f states for the two defects, with the results for the SVBT V defect shown in Fig. 9(c). Also here, we note that the electronic properties are in qualitative agreement for the two methods. However, as seen in Fig. 9(c), the PBE+U predicts a splitting of the two occupied Ce4f states for the SVBT V defect, which does not appear in the DOS generated from the DFTB+U calculation. Overall, the DFTB+U method reproduces the properties of the oxygen vacancy at the (110) surface well compared to our reference method.

22


Page 22 of 55

Page 23 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Surface steps at CeO2 (111) – 1D Surface steps are common one dimensional defects and can constitute a large fraction of active sites on a catalytic surface. It is not our intention here to fully investigate all possible structures of steps at the ceria surface. Instead, we consider one step structure as an example, namely the so-called T ype I step, which has previously been found to constitute the most stable step type in the DFT study of Kozlov et al. 60 The step model is illustrated in Fig. 10. The step formation energy, Estep , is the energy penalty to create the step from the relaxed extended (111) surface. To calculate Estep , we used the same procedure as in Ref. 60 in which the step energy is obtained by fitting the following energy expression to a set of calculated total energies for stepped slabs of different thicknesses and terrace lengths (different step–step separations):

Estepped slab = N · Ebulk + 2wL · Esurf + 2w · Estep .

(8)

In Eq. 8, N is the number of formula units in the slab supercell, Ebulk is the energy of a bulk formula unit, L is the length of the slab perpendicular to the step, Esurf is the surface energy, Estep is the step formation energy per unit length of the step, and w is the width of the slab along the step. Slabs consisting of 3, 4 and 5 triple-layers, and with terrace lengths of 23.3, 31.0 and 38.8 ˚ A were used. The fitted step energies Estep become 0.84 and 1.25 eV/nm, respectively, using DFTB+U and PBE+U . Thus, here the absolute agreement is similar to that of the surface energy for the (111) surface as calculated in section , but the DFTB+U now underestimates the formation energy. However, as in the case of the (111) surface, we will demonstrate in the following that the local structure, i.e. the relaxation, and the vacancy formation energies at the T ype I step are in better agreement. To quantify the relaxation imposed by the step, we have measured some key geometrical parameters in the optimised DFTB+U and PBE+U structures. These parameters are

23



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

defined in Fig. 10(a) and compared in Table 6. The agreement is very good. We have also calculated the vacancy formation energy at the step by removing an oxygen ion close to the lower terrace as was proposed in the DFT study of Kozlov et al. 61 (cf. Fig. 10(a)). The vacancy formation energy is given in Table 2. Here, DFTB+U overestimates the formation energy compared to the PBE+U results. We furthermore analysed the electron DOS and the DOS projected on the Ce3+ ions, cf. Fig. 10(b). The overall agreement with the two methods is good. However there is a slight discrepancy when it comes to the position of the localised Ce4f state.

Nano structures – 0D In this section, we test the performance of the generated DFTB+U parameters when applied to ceria nano structures. We consider stoichiometric nanoparticles (NPs) as well as two types of partially reduced NPs, namely inherently under-stoichiometric NPs without explicit oxygen vacancies and stoichiometric particles with explicit oxygen vacancies. 11 Stoichiometric and under-stoichiometric ceria nanoparticles In experiments, the most frequently occurring shape for ceria NPs is the octahedral one, exposing almost exclusively the most stable (111) facet, c.f. Fig.11(a). This is in agreement with theoretical modelling using the Wulff construction scheme. Tetrahedrally shaped NPs have not been observed experimentally but have been suggested, based on DFT calculations, to constitute the most stable shape for particles smaller than 2 nm in diameter. 62 Besides the octahedrally shaped NPs, ceria cubes have also been synthesized. 5 All simulations of nanoparticles were simulated suing large supercells spanning from 25 × 25 × 25 ˚ A for the smallest particles to 35 × 35 × 35 ˚ A for the largest ones. All strucutres were fully optimized without any constraints using both VASP and DFTB+. Perfect octahedrally shaped CeO2 NPs (Oct-NPs) cut from the fluorite bulk structure are always cerium rich and consequently contain a number of Ce3+ ions. Previous DFT 24


Page 24 of 55

Page 25 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


studies suggested that these Ce3+ -ions are evenly distributed over the Oct-NP edges, 63 cf. Fig. 11(a). The number of Ce3+ -ions in the particle is proportional to the side length, i.e the number of Ce3+ grows linearly with particle diameter. As suggested by Migani et al., 64 stoichiometric particles can be constructed by truncation of the NP edges, i.e. removal of Ce-ions exposing (100) planes, see Fig.11(b). The number of truncating planes and their individual sizes depend on the size of the particle and can be obtained by analysis of the number series governing the number of Ce and O ions in octahedral particles. Perfect tetrahedral nano particles exposing only (111) facets are always oxygen rich. In the current investigation we use a similar procedure as that leading to the stoichiometric Oct-NPs in which stoichiometric tetrahedra are obtained by truncating some or all of the corners with (111) planes, cf. Fig.11(c). Just as for the octahedral particles, the number of truncating planes and their sizes depend on the particle size and can be obtained from the number series governing the number of Ce and O in tetrahedral particles. Perfect cubes, exposing only (100) facets, cut out from the fluorite structure exhibit either oxygen rich or cerium rich faces. Similar to the (100) surface described earlier in this study, we have used stoichiometric particles that were obtained by introducing oxygen vacancies at the surface of completely oxygen terminated cubes, Fig. 11(d). The positions of the vacancies were determined from a point-charge model using a Monte-Carlo (MC) procedure. In our MC approach, vacancies where allowed to hop on lattice sites on the surface of the cubes and the total electrostatic energy was calculated for each step and used to determine whether a hop should be allowed or not. For each particle, at least 40000 MC steps where used. This procedure led to particles with vacancies spread evenly over the surface and the vacancy pattern resembles the one used for the reconstructed (100)r surface presented in Sec. (d). For each type of NPs, both stoichiometric and reduced, we have studied NPs of different sizes and calculated their stabilities compared to bulk CeO2 in the fluorite structure. Fig. 12 display the results for particles with up to 200 formula units, calculated with both

25



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

DFTB+U and PBE+U . The stabilities calculated with the PBE+U method are clearly very well reproduced by our DFTB+U method, both between the different shapes (octahedral, tetrahedral and cube) and the different stoichiometries (perfect octahedral and truncated octahedral). Oxygen vacancies in truncated octahedrally shaped ceria nano particles Vacancy formation energies in ceria nanoparticles has been shown by DFT calculations to be significantly lower than those of the extended surfaces. 64 Here, we do not explore every possibile location of the oxygen vacancy on the ceria NP for all sizes. Instead, we compare the formation energy, Evac , for a vacancy formed at the Ce40 O80 cluster surface using the structure presented in Ref. 64 in which the vacancy is located at one of the small (100) surface facets, see Fig. 13(a). With our PBE+U approach, we obtain an Evac value of 0.60 eV. This value is slightly lower than that presented in Ref. 64 which used a similar DFT approach albeit one with a lower U -value (4 eV), and a different GGA functional (PW91). With DFTB+U we obtain an Evac value of 0.47 eV in close agreement with our PBE+U reference. The electronic DOS, obtained with the two methods, are compared in Fig. 13(b). The agreement is, as for the other cases tested, very good and the downshift of the occupied Ce4f state when located at a corner site compared to e.g. at the (111)-surface is well captured by the DFTB+U method.

Summary of validation tests We now summarise the performance of the generated parameters for the calculated properties of stoichiometric and partially reduced ceria. Overall, considering all calculated properties and all the ceria systems considered here, the DFTB+U parametrisation presented here demonstrate we good qualitative agreement with our reference PBE+U method and often good quantitative agreement as well. It is important to note that some of the larger errors are a consequence of a deliberate sacrifice 26


Page 26 of 55

Page 27 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


in order to improve the descriptions of other properties. One important example is that of the relative stabilities of the rutile and the fluorite phases of CeO2 , where we intentionally decreased the stability of the rutile phase (but not too much) to achieve a better description of the fluorite (111) surface. The structural discrepancies between the two methods are small. Here, the average error with respect to measured bond lengths is small, ±0.4 % (excluding the rutile structure), with the largest one corresponding to 0.05 ˚ A. For relative relaxation effects, the percentage error varies more, but also here the maximum error is small, being on the order of 0.05 ˚ A. For all considered reduction energies, our new parameter set gives an average error of ±11%. All calculated vacancy formation energies reported in Table 2 are displayed in Fig.14. In the figure each box represents a system type, the qualitative agreement between the two methods is illustrated by the good correspondence between similar boxes for the two methods. We conclude that we are able to describe vacancy formation energies in ceria systems of various dimensions. Individual DFTB+U agree with the PBE+U reference values within, typically, 0.1–0.4 eV. The largest discrepancy is found for the T ype I step and the SVBT V vacancies. We note that these form very similar local structures after the vacancy formation, namely a twocoordinated O atom. Regarding the accuracy when it comes to absolute numbers, we are somewhat limited by the pair-wise repulsive potential we use, which might be too simplistic an approach to enable broad transferability between situations involving different ionic coordinations. For example, improving on the bulk vacancy formation energy has in our attempts led to a worsening of those of the surface ones, and vice versa. Thus, within the current formalism of SCC-DFTB, we are satisfied with the qualitative agreement, and importantly, with reduction energies that are of the right magnitude when compared to the reference method. For the electronic properties, the error with respect to our reference method varies for the different partially reduced systems with errors as large as 0.4 eV. However, it is important to

27



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

note that qualitatively, the trends for the different systems persists for all studied properties. Furthermore, the parameters also give qualitative structural agreement when it comes to the electron localisation, i.e. Ce3+ formation in the different cases, which has been one of the goals with the parameterisation.

Conclusion We have presented DFTB+U parameters for Ce–O interactions with Ce in both the +III and +IV formal oxidation states, which are applicable to studies of stoichiometric ceria (CeO2 ), partially reduced ceria (CeO2−x ) and fully reduced Ce2 O3 . The new parameters have been thoroughly validated for bulk, surface and nanoparticle ceria systems of different stoichiometry by comparison with PBE+U calculations, which constitutes reference method which we want our new method to reproduce. The validation tests showed that we are indeed able to describe reduction reactions in ceria of different dimensionality. For all studied properties we obtain a qualitative, and often quantitative, agreement with our reference DFT calculations. The targeted tuning of the parameters unfortunately leads to reduced transferability, where, for example, the Ce metal is not well described. We also present guidelines for an iterative SCC-DFTB parameterisation scheme for reducible oxides, which was formulated such that particular emphasis was placed on attaining good to improve the accuracy (i.e. agreement with PBE+U ) for calculated reduction energies. In particular, we increased the number of degrees of freedom in the parameterisation by optimising parameters in the SCC-DFTB formalism that are usually obtained from atomic reference calculations. In conclusion, the generated parameters were shown to provide a robust, accurate, and computationally efficient alternative in theoretical studies of ceria systems with the possibility of reaching (at least) an order of magnitude larger size and time scales compared to DFT

28


Page 28 of 55

Page 29 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


simulations. Further development with additional SCC-DFTB parameters for Ce interacting with other atom types will open up new possiblities for ceria systems.

Acknowledgements This work was supported by the the Swedish Research Council (VR), the ˚ Angpanneföreningens forskningsstiftelse (˚ Aforsk), and STINT (The Swedish Foundation for International Cooperation in Research and Higher Education). Funding from the National Strategic e-Science program eSSENCE is greatly acknowledged. This work was also supported by the COST Action CM1104 “Reducible oxide chemistry, structure and functions”. The simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at UPPMAX and NSC. JK would like to acknowledge the Hanse-Wissenschaftskolleg (HWK) for support.

Supporting Material Content Supporting material compare the quallity of the parameter set presented in the main text to an alternative parameter set where the rutile phase is better described. The key properties compared are: • Relative stability of CeO2 in the flourite and rutile structure. • Surface energies for the low index surfaces of CeO2 . • Vacancy formation energies in the bulk and on the low index surfaces of CeO2 .

References (1) Trovarelli, A. In Catalysis by Ceria and Related Materials; Hutchings, G. J., Ed.; Imperial College Press, 2002.

29



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(2) Park, S.; Vohs, J. M.; Gorte, R. J. Direct Oxidation of Hydrocarbons in a Solid-Oxide Fuel Cell. Nature 2000, 404, 265–267. (3) Paier, J.; Penschke, C.; Sauer, J. Oxygen Defects and Surface Chemistry of Ceria: Quantum Chemical Studies Compared to Experiment. Chem. Rev. 2013, 113, 3949– 3985. (4) Zhang, F.; Jin, Q.; Chan, S.-W. Ceria Nanoparticles: Size, Size Distribution, and Shape. J. Appl. Phys. 2004, 95, 4319–4326. (5) Mai, H.-X.; Sun, L.-D.; Zhang, Y.-W.; Si, R.; Feng, W.; Zhang, H.-P.; Liu, H.-C.; Yan, C.-H. Shape-Selective Synthesis and Oxygen Storage Behavior of Ceria Nanopolyhedra, Nanorods, and Nanocubes. J. Phys. Chem. B 2005, 109, 24380–24385. (6) Zhou, K.; Wang, X.; Sun, X.; Peng, Q.; Li, Y. Enhanced Catalytic Activity of Ceria Nanorods From Well-Defined Reactive Crystal Planes. J. Catal. 2005, 229, 206 – 212. (7) Zhang, J.; Ohara, S.; Umetsu, M.; Naka, T.; Hatakeyama, Y.; Adschiri, T. Colloidal Ceria Nanocrystals: A Tailor-Made Crystal Morphology in Supercritical Water. J. Adv. Mater. 2007, 19, 203–206. (8) Yang, Z.; Zhou, K.; Liu, X.; Tian, Q.; Lu, D.; Yang, S. Single-crystalline Ceria Nanocubes: Size-Controlled Synthesis, Characterization and Redox Property. Nanotechnology 2007, 18, 185606. (9) Wu, Z.; Li, M.; Howe, J.; Meyer, H. M.; Overbury, S. H. Probing Defect Sites on CeO2 Nanocrystals with Well-Defined Surface Planes by Raman Spectroscopy and O2 Adsorption†. Langmuir 2010, 26, 16595–16606. (10) Xu, J.; Harmer, J.; Li, G.; Chapman, T.; Collier, P.; Longworth, S.; Tsang, S. C. Size Dependent Oxygen Buffering Capacity of Ceria Nanocrystals. Chem. Commun. 2010, 46, 1887–1889. 30


Page 30 of 55

Page 31 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(11) Kullgren, J.; Hermansson, K.; Broqvist, P. Supercharged Low-Temperature Oxygen Storage Capacity of Ceria at the Nanoscale. J. Phys. Chem. Lett. 2013, 4, 604–608. (12) Renuka, N. K.; Harsha, N.; Divya, T. Supercharged Ceria Quantum Dots With Exceptionally High Oxygen Buffer Action. RSC Adv. 2015, 5, 38837–38841. (13) Silva, G. A. Nanomedicine: Seeing the Benefits of Ceria. Nat. nanotech. 2006, 1, 92–94. (14) Chen, J.; Patil, S.; Seal, S.; McGinnis, J. F. Rare Earth Nanoparticles Prevent Retinal Degeneration Induced by Intracellular Peroxides. Nature Nanotechnology 2006, 1, 142–150. (15) Balducci, G.; Kaˇspar, J.; Fornasiero, P.; Graziani, M.; Islam, M. S. Surface and Reduction Energetics of the CeO2-ZrO2 Catalysts. J. Phys. Chem. B 1998, 102, 557–561. (16) Burbano, M.; Marrocchelli, D.; Yildiz, B.; Tuller, H. L.; Norberg, S. T.; Hull, S.; Madden, P. A.; Watson, G. W. A Dipole Polarizable Potential for Reduced and Doped CeO2 Obtained From First Principles. J. Phys.: Condens. Matter 2011, 23, 255402. (17) Conesa, J. Computer Modeling of Surfaces and Defects on Cerium Dioxide. Surf. Sci. 1995, 339, 337–352. (18) Gotte, A.; Sp˚ angberg, D.; Hermansson, K.; Baudin, M. Molecular Dynamics Study of Oxygen Self-Diffusion in Reduced CeO2. Solid State Ionics 2007, 178, 1421–1427. (19) Gotte, A.; Hermansson, K.; Baudin, M. Molecular Dynamics Simulations of Reduced CeO2: Bulk and Surfaces. Surf. Sci. 2004, 552, 273–280. (20) Sayle, T. X. T.; Parker, S. C.; Catlow, C. R. A. The Role of Oxygen Vacancies on Ceria Surfaces in the Oxidation of Carbon Monoxide. Surf. Sci. 1994, 316, 329–336. (21) Vyas, S.; Grimes, R. W.; Gay, D. H.; Rohl, A. L. Structure, Stability and Morphology of Stoichiometric Ceria Crystallites. J. Chem. Soc., Faraday Trans 1998, 94, 427–434. 31



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(22) van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard, W. A. ReaxFF: A Reactive Force Field for Hydrocarbons. J. Phys. Chem. A 2001, 105, 9396–9409. (23) van Duin, A. C. T.; Strachan, A.; Stewman, S.; Zhang, Q.; Xu, X.; Goddard, W. A. ReaxFFSiO Reactive Force Field for Silicon and Silicon Oxide Systems. J. Phys. Chem. A 2003, 107, 3803–3811. (24) Zhang, Q.; C ¸ aˇgın, T.; van Duin, A.; Goddard, W. A.; Qi, Y.; Hector, L. G. Adhesion and Nonwetting-Wetting Transition in the Al/α - Al2O3 Interface. Phys. Rev. B 2004, 69, 045423. (25) Broqvist, P.; Kullgren, J.; Wolf, M. J.; van Duin, A. C. T.; Hermansson, K. ReaxFF Force-Field for Ceria Bulk, Surfaces, and Nanoparticles. The J. Phys. Chem. C 2015, 119, 13598–13609. (26) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-Consistent-Charge Density-Functional Tight-Binding Method for Simulations of Complex Materials Properties. Phys. Rev. B 1998, 58, 7260–7268. (27) Frauenheim, T.; Seifert, G.; Elsterner, M.; Hajnal, Z.; Jungnickel, G.; Porezag, D.; Suhai, S.; Scholz, R. A Self-Consistent Charge Density-Functional Based Tight-Binding Method for Predictive Materials Simulations in Physics, Chemistry and Biology. physica status solidi (b) 2000, 217, 41–62. (28) Moreira, N. H.; Dolgonos, G.; Aradi, B.; da Rosa, A. L.; Frauenheim, T. Toward an Accurate Density-Functional Tight-Binding Description of Zinc-Containing Compounds. J. Chem. Theory Comput. 2009, 5, 605–614. (29) Dolgonos, G.; Aradi, B.; Moreira, N. H.; Frauenheim, T. An Improved Self-ConsistentCharge Density-Functional Tight-Binding (SCC-DFTB) Set of Parameters for Simulation of Bulk and Molecular Systems Involving Titanium. J. Chem. Theory Comput. 2010, 6, 266–278. 32


Page 32 of 55

Page 33 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(30) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865–3868. (31) Blöchl, P. E. Projector Augmented-Wave Method. Phys. Rev. B 1994, 50, 17953–17979. (32) Dudarev, S. L.; Botton, G. A.; Savrasov, S. Y.; Humphreys, C. J.; Sutton, A. P. Electron-Energy-Loss Spectra and the Structural Stability of Nickel Oxide:

An

LSDA+U Study. Phys. Rev. B 1998, 57, 1505–1509. (33) Andersson, D. A.; Simak, S. I.; Johansson, B.; Abrikosov, I. A.; Skorodumova, N. V. Modeling of CeO2, Ce2O3, and CeO2-x in the LDA+U Formalism. Phys. Rev. B 2007, 75, 035109. (34) Castleton, C. W. M.; Kullgren, J.; Hermansson, K. Tuning LDA + U for Electron Localization and Structure at Oxygen Vacancies in ceria. J. Chem. Phys 2007, 127, 244704. (35) Da Silva, J. L. F.; Ganduglia-Pirovano, M. V.; Sauer, J.; Bayer, V.; Kresse, G. Hybrid Functionals Applied to Rare-Earth Oxides: The Example of Ceria. Phys. Rev. B 2007, 75, 045121. (36) Loschen, C.; Carrasco, J.; Neyman, K. M.; Illas, F. First-Principles LDA+U and GGA+U Study of Cerium Oxides: Dependence on the Effective U Parameter. Phys. Rev. B 2007, 75, 035115. (37) Jiang, Y.; Adams, J. B.; Schilfgaarde, M. v. Density-functional Calculation of CeO2 Surfaces and Prediction of Effects of Oxygen Partial Pressure and Temperature on Stabilities. J. Chem. Phys 2005, 123, 064701. (38) Kresse, G.; Hafner, J. Ab Initio Molecular Dynamics for Liquid Metals. Phys. Rev. B 1993, 47, 558–561.

33



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(39) Kresse, G.; Hafner, J. Ab Initio Molecular-Dynamics Simulation of the Liquid-Metal– Amorphous-Semiconductor Transition in Germanium. Phys. Rev. B 1994, 49, 14251– 14269. (40) Kresse, G.; Furthm¨ uller, J. Efficiency of Ab-Initio Total Energy Calculations for Metals and Semiconductors Using a Plane-Wave Basis Set. Computational Materials Science 1996, 6, 15–50. (41) Kresse, G.; Furthm¨ uller, J. Efficient Iterative Schemes for Ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set. Phys. Rev. B 1996, 54, 11169–11186. (42) Hourahine, B.; Sanna, S.; Aradi, B.; Köhler, C.; Niehaus, T.; Frauenheim, T. SelfInteraction and Strong Correlation in DFTB†. J. Phys. Chem. A 2007, 111, 5671–5677. (43) Aradi, B.; Hourahine, B.; Frauenheim, T. DFTB+, a Sparse Matrix-Based Implementation of the DFTB Method†. J. Phys. Chem. A 2007, 111, 5678–5684. (44) Bodrog, Z.; Aradi, B.; Frauenheim, T. Automated Repulsive Parametrization for the DFTB Method. J. Chem. Theory Comput. 2011, 7, 2654–2664. (45) Wahiduzzaman, M.; Oliveira, A. F.; Philipsen, P.; Zhechkov, L.; van Lenthe, E.; Witek, H. A.; Heine, T. DFTB Parameters for the Periodic Table: Part 1, Electronic Structure. J. Chem. Theory Comput. 2013, 9, 4006–4017. (46) Chou, C.-P.; Nishimura, Y.; Fan, C.-C.; Mazur, G.; Irle, S.; Witek, H. A. Automatized Parameterization of DFTB Using Particle Swarm Optimization. J. Chem. Theory Comput. 2016, 12, 53–64. (47) Hellström, M.; Jorner, K.; Bryngelsson, M.; Huber, S. E.; Kullgren, J.; Frauenheim, T.; Broqvist, P. An SCC-DFTB Repulsive Potential for Various ZnO Polymorphs and the ZnO–Water System. The J. Phys. Chem. C 2013,

34


Page 34 of 55

Page 35 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(48) Gale, J. D.; Rohl, A. L. The General Utility Lattice Program (GULP). Mol. Simul. 2003, 29, 291–341. (49) Liu, W.; Dolg, M. Benchmark Calculations for Lanthanide Atoms: Calibration of Ab Initio and Density-Functional Methods. Phys. Rev. A 1998, 57, 1721–1728. (50) Ren, C.-Y. Relativistic Density-Functional All-Electron Calculations of Interconfigurational Energies of Lanthanide Atoms. J. Chem. Phys 2004, 121, 11073–11082. (51) Jones, R. O.; Gunnarsson, O. The Density Functional Formalism, its Applications and Prospects. Rev. Mod. Phys. 1989, 61, 689–746. (52) Harris, J.; Jones, R. O. Density Functional Theory of 3d-Transition Element Atoms. J. Chem. Phys 1978, 68, 3316–3317. (53) Furche, F.; Perdew, J. P. The Performance of Semilocal and Hybrid Density Functionals in 3d Transition-Metal Chemistry. J. Chem. Phys 2006, 124, 044103. (54) Averill, F. W.; Painter, G. S. Steepest-Descent Determination of Occupation Numbers and Energy Minimization in the Local-Density Approximation. Phys. Rev. B 1992, 46, 2498–2502. (55) Bärnighausen, H.; Schiller, G. The Crystal Structure of A-Ce2O3. Journal of the Less Common Metals 1985, 110, 385–390. (56) Morss, L. R. Comparative Thermochemical and Oxidation-Reduction Properties of Lanthanides and Actinides. Handbook on the Physics and Chemistry of Rare Earths 1994, 18, 239–291. (57) Nolan, M.; Grigoleit, S.; Sayle, D. C.; Parker, S. C.; Watson, G. W. Density Functional Theory Studies of the Structure and Electronic Structure of Pure and Defective Low Index Surfaces of Ceria. Surf. Sci. 2005, 576, 217–229.

35



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(58) Galea, N. M.; Scanlon, D. O.; Morgan, B. J.; Watson, G. W. A GGA+U Study of the Reduction of Ceria Surfaces and Their Partial Reoxidation Through NO2 Adsorption. Mol. Simul. 2009, 35, 577–583. (59) Kullgren, J.; Hermansson, K.; Castleton, C. Many Competing Ceria (110) Oxygen Vacancy Structures: From Small to Large Supercells. J. Chem. Phys 2012, 137, 044705. (60) Kozlov, S. M.; Vi˜ nes, F.; Nilius, N.; Shaikhutdinov, S.; Neyman, K. M. Absolute Surface Step Energies: Accurate Theoretical Methods Applied to Ceria Nanoislands. J. Phys. Chem. Lett. 2012, 3, 1956–1961. (61) Kozlov, S. M.; Neyman, K. M. O Vacancies on Steps on CeO2(111) Surface. Phys. Chem. Chem. Phys. 2014, 16, 7823–7829. (62) Migani, A.; Neyman, K. M.; Bromley, S. T. Octahedrality Versus Tetrahedrality in Stoichiometric Ceria Nanoparticles. Chem. Commun. 2012, 48, 4199–4201. (63) Loschen, C.; Migani, A.; Bromley, S. T.; Illas, F.; Neyman, K. M. Density functional studies of model cerium oxide nanoparticles. Phys. Chem. Chem. Phys. 2008, 10, 5730– 5738. (64) Migani, A.; Vayssilov, G. N.; Bromley, S. T.; Illas, F.; Neyman, K. M. Greatly Facilitated Oxygen Vacancy Formation in Ceria Nanocrystallites. Chem. Commun. 2010, 46, 5936–5938.

36


Page 36 of 55

Page 37 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1: Calculated cellparameters, bulk modulus and electronic energy gaps for bulk fluorite CeO2 using DFTB+U and PBE+U . The value in parenthesis refers to the energy gap calculated from the O2p valence band edge to the center of the unoccupied Ce4f band. a (˚ A) B (GPa)

DFTB+U 5.48 182

PBE+U 5.49 179

∆ǫO2p→Ce4f (unocc) (eV) ∆ǫO2p→Ce5d(unocc) (eV)

2.58 (2.95) 5.51

2.42 (2.86) 5.37

37



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 2: Oxygen vacancy formation energies, Evac , in eV for various ceria structures calculated with DFTB+U and PBE+U according to formula 2. The reduction energy for Ce2 O3 according to formula 4. This energy is given in eV per half O2 product to comply with the definition of the vacancy formation energies. System Bulk VN N Bulk VN N N

DFTB+U 2.93 3.04

PBE+U 3.26 3.22

CeO2 → A-Ce2 O3

3.18

2.30

(111) (111) (111) (111)

2.38 2.38 2.39 2.21

2.58 2.34 2.43 2.12

(110) SVIP V (110) SVBT V

1.97 1.82

1.94 1.56

(111) T ype 1 step

1.91

1.48

(100)@Ce40 O80

0.47

0.60

SVN N SVN N N SSVN N SSVN N N

38


Page 38 of 55

Page 39 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 3: Structural relaxation around an oxygen vacancy in fluorite bulk CeO2 and position of occupied Ce4f state relative to the O2p band maximum. Distances are in ˚ A and energy difference in eV. The atomic labels are decoded in Fig. 5(a) and 5(b). VN N rCe1 −Ce2 rCe1 −Ce3 rCe1 −Ce4 rCe2 −Ce3 rCe2 −Ce4 rCe3 −Ce4 O2p→Ce4f (occ) ∆ǫ VN N N rCe1 −Ce2 rCe1 −Ce3 rCe1 −Ce4 rCe2 −Ce3 rCe2 −Ce4 rCe3 −Ce4 O2p→Ce4f (occ) ∆ǫ

DFTB+U

PBE+U

4.18 4.18 4.18 4.12 4.12 4.05 1.99

4.20 4.13 4.13 4.13 4.13 4.07 1.52

4.13 4.15 4.15 4.13 4.13 4.14 2.02

4.12 4.15 4.15 4.13 4.13 4.11 1.77

39



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 4: Calculated surface energies (J/m2 ), surface relaxation ρ and surface rumpling ǫ in % (see text for details) using DFTB+U and PBE+U . The (111), (110) and (110)r surfaces are schematically illustrated in Fig. 7. DFTB+U CeO2 (111) Esurf 0.93 ρ -0.7 ǫ -0.8

PBE+U 0.71 -0.9 -0.9

CeO2 (110) Esurf 1.16 ρ -7.9 ǫ +18.5

1.09 -7.0 +13.0

CeO2 (100)r Esurf 1.57 ρ -6.1 ǫ -14.0

1.54 -5.0 -14.0

40


Page 40 of 55

Page 41 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 5: Structural deformation around an oxygen vacancy on the CeO2 (111) surface. Distances are in ˚ A. Notations are explained in Fig. 8 and in the text. DFTB+U

PBE+U

3.79 3.95 4.09 3.94

3.83 3.95 4.14 3.96

3.59 4.23 3.65 4.01

3.61 4.25 3.67 4.01

0.21

0.16

0.30

0.29

SVN N rO1 −O2 rO2 −O3 rO3 −O4 rO4 −O5 SVN N N rO1 −O2 rO2 −O3 rO3 −O4 rO4 −O5 SSVN N ∆z(ON N – ON N N ) SSVN N N ∆z(ON N – ON N N )

41



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 6: Step formation energies and relaxation at a type I step on the CeO2 (111) surface calculated with DFTB+U and PBE+U .

d1 d2 d3 θ1 θ2

(˚ A) (˚ A) (˚ A) (◦ ) (◦ )

Estep (eV/nm)

DFTB+U 3.23 5.78 3.23 148.4 148.8

PBE+U 3.24 5.74 3.28 147.4 147.3

0.84

1.25

42


Page 42 of 55

Page 43 of 55


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment



Page 44 of 55

Page 45 of 55





Page 46 of 55

Page 47 of 55





Page 48 of 55

Page 49 of 55





Page 50 of 55

Page 51 of 55





Page 52 of 55

Page 53 of 55





Page 54 of 55

Page 55 of 55



Self-Consistent-Charge Density-Functional Tight-Binding (SCC-DFTB

Recommend Documents