DFT-based Monte Carlo Simulations of Impurity ... - ACS Publications

Article pubs.acs.org/JPCC

DFT-based Monte Carlo Simulations of Impurity Clustering at CeO2(111) Jolla Kullgren,*,† Matthew J. Wolf,† Pavlin D. Mitev,† Kersti Hermansson,† and Wim J. Briels‡,§ †

Department of Chemistry − Ångström Laboratory, Uppsala University, Box 538, 751 21 Uppsala, Sweden Computational Biophysics, University of Twente, Box 217, 7500 AE Enschede, The Netherlands § Forschungszentrum Jülich, ICS 3, D-52425 Jülich, Germany ‡

S Supporting Information *

ABSTRACT: The interplay between energetics and entropy in determining defect distributions at ceria(111) is studied using a combination of DFT+U and lattice Monte Carlo simulations. Our main example is fluorine impurities, although we also present preliminary results for surface hydroxyl groups. A simple classical force-field model was constructed from a training set of DFT+U data for all symmetrically inequivalent (F−)n(Ce3+)n nearestneighbor clusters with n = 2 or 3. Our fitted model reproduces the DFT energies well. We find that for an impurity concentration of 15% at 600 K, straight and hooked linear fluorine clusters are surprisingly abundant, with similarities to experimental STM images from the literature. We also find that with increasing temperature the fluorine cluster sizes show a transition from being governed by an attractive potential to being governed by a repulsive potential as a consequence of the increasing importance of the entropy of the Ce3+ ions. The distributions of surface hydroxyl groups are noticeably different.

■

investigations6,8 which found surface oxygen vacancy dimers to be unstable relative to well separated single vacancies. Given the fact that none of the current state-of-the-art DFT methods for ceria predicts clustering of surface oxygen vacancies, several investigators have searched for alternative explanations. In a 2014 paper,9 inspired by the detection of concentrations of fluorine impurities of up to 10% in ceria samples from the same source as those used in the STM experiments,1−4 we proposed that the observed defects may, in fact, have been fluorine impurities. Our DFT calculations at the PBE+U level showed that fluorine dimers in the top (111) surface layer are more stable than two isolated fluorine impurities by ∼0.04 eV per fluorine. A later DFT study also showed that fluorine should not readily desorb.7 Moreover, we demonstrated9 that the STM images of isolated fluorine, hydroxyl, and oxygen vacancy defects are essentially identical. Following that paper, Wu and Gong10 carried out DFT calculations of surface fluorine trimer clusters and found that the various linear, bent, and triangular arrangements tested were all close in energy to each other and close to that of an “F dimer + F monomer”, with a small preference for a compact, triangular structure. The authors concluded that “linear configurations of fluorine impurities will not dominate the

INTRODUCTION

The physical and chemical properties of technologically relevant oxide materials often depend critically on the nature, distribution, and properties of defects. Thus defect identification and characterization is the subject of much research effort. Analysis based solely on the results of experimental investigations is often not sufficient, and theoretical calculations play an increasingly important role in supporting and explaining experimental findings. The present study discusses the identity and properties of point defects on the most stable surface of ceria, namely, CeO2(111). There is an ongoing discussion in the literature concerning the tendency (or lack thereof) for oxygen vacancies to form stable clusters at this surface. Numerous experimental investigations of CeO2(111) using scanning tunneling microscopy (STM) have been undertaken for single-crystal samples.1−4 In these studies, the surface defects were assumed to be oxygen vacancies and were observed to form nearestneighbor clusters of various shapes, such as triangular and extended linear clusters, with an apparent preference for the latter, at least in ref 4. Interestingly, very few theoretical investigations supporting these observations and interpretations have been presented. Early classical force-field studies by Conesa5 suggested that surface oxygen vacancies would associate into pairs on the ceria(111) surface. However, these results were later refuted by density functional theory (DFT) © XXXX American Chemical Society

Received: January 11, 2017 Revised: April 18, 2017 Published: May 2, 2017 A

DOI: 10.1021/acs.jpcc.7b00299 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 1. Top view of the CeO2 (111) surface, with the fluorine defect clusters used in the generation of the DFT-based force field. The Ce3+ ions are shown at their lowest-energy locations (up to symmetric equivalence).

surface, which is largely conflicted with the STM observations”. As an alternative explanation they suggested that subsurface hydrogen impurities assist in the formation of linear surface oxygen vacancy trimers. In our opinion, the relevance of the Hassisted mechanism suggested by Wu and Gong is questionable as it turns out that their clusters are in fact highly unstable with respect to decomposition into isolated defects.11 It is striking that in the arguments used in the literature so far for ceria surface defects only energetic considerations were taken into account. Because energy differences involved in the clustering correspond to temperatures of a few hundred to 1000 K, it is clear that entropic effects should play an important role at the temperatures of the experiments. Moreover, the Ce3+ ions associated with the defects are often placed at their most stable positions in calculations, but because the energy required to bring the Ce3+ ions (or rather the electrons) to less favorable positions again is on the order of several hundreds of K, their contributions to the total entropy should also be important. Here we set forth to investigate, by means of Monte Carlo (MC) simulations, the influence of entropy on the type of ceria surface defect structures discussed above. The MC simulations are based on an interaction model derived from DFT. Similar approaches have been successfully used to study dopant distributions in bulk ceria with oxygen vacancies.12,13 We restrict ourselves to 2D ceria(111) systems containing surface fluorine impurities and the corresponding Ce3+ ions in the layer underneath. The restriction to a 2D system is motivated by the fact that in fluorine-doped CeO2 it costs 0.48 eV9 to bring a fluorine impurity from the surface to the first subsurface layer. We will demonstrate that this system and our simulation approach reveals some interesting physics; in particular, we will find that at higher temperatures linear structures become quite abundant, even though compact structures are favored by energetic considerations alone. At the end of the paper we will briefly discuss the relevance of our results to the interpretation of the STM measurements mentioned above.

Figure 2. Correlation plot between the formation energies (per F atom) of the 78 clusters in the training set, calculated with the two different methods. The DFT energies from the training set are the values on the x axis and the energies of the fitted model function (from expressions 1, 2, and 5 in the text) are given on the y axis. For the reader’s convenience, two energy scales are given in the Figure.

several simulations using p(30 × 30) and p(20 × 20) supercells and consistently found the results to be robust. The simulations were done with the GULP code,14 in such a way that the trial moves consisted of swaps between any species on a randomly chosen sublattice, that is, the cation or the anion lattice. Acceptance probabilities were calculated according to the usual Metropolis criterion based on the Boltzmann factor associated with the energy difference between the trial and current state. One should note that the fluoride anions and Ce3+ cations reside on different sublattices and that no hops between them were allowed and no double occupancy of sites was allowed. In each Monte Carlo simulation the energy was monitored. Properties were sampled from one hundred thousand accepted moves after the energy had become constant apart from random fluctuations. To check for unnoticeable slow condensation processes, several of the runs were started in two ways, namely, once from a random configuration and once from a condensed configuration; in all cases the system developed toward the same final state. Interaction Model. A simple classical interaction model was used of the following mathematical form

■

METHODS Monte Carlo Setup. A top-view of CeO2(111) is shown in Figure 1. The structural model in the Monte Carlo simulations consisted of the top two atomic layers of the (111) facet, as mentioned in the Introduction. In the top layer, a set number of oxide ions (white circles) were replaced by fluorine ions (green circles), while in the cation layer immediately below it an equal number of Ce4+ ions (black circles) were changed to Ce3+ ions (red circles). The 3+/4+ charges attributed to the cerium ions are nominal charges serving only to distinguish the two types of cerium ions. All simulations were carried out under 2D periodic boundary conditions, using a p(40 × 40) surface supercell. We repeated

Φmodel =

∑

ϕ(a , b) (1)

where ∑ is the sum over all unique pairs of fluorines or Ce3+ in the system, and the pair contributions are given by B


Figure 3. Monte Carlo results at three different temperatures (rows 2-4 from the top) and with three different fluorine defect concentrations. In each case, the left panel shows the average fraction of F− ions taking part in clusters ranging from monomers to hexamers. The right panel shows a snapshot (with the F− ions marked in black) chosen such that its distribution of F participation in clusters of different sizes is very similar to the run-time average shown to the left of it. With each cluster size the contribution of open clusters, as defined in the main text, is indicated by white boxes. The panels in the top row show similar histograms and snapshots obtained by distributing the F− ions randomly.

The Journal of Physical Chemistry C Article

C


Article


figure displays the Ce3+ ions in their lowest energy positions, but in fact all possible unique localization patterns of the associated Ce3+ ions were considered, with the restriction that each Ce3+ ion be located on a site which is a nearest-neighbor to at least one of the fluorine ions. Altogether our training set consisted of 78 DFT+U defect cluster formation energy values, which are given in the Supporting Information. All DFT calculations were performed with the Vienna Ab-initio Simulation Package (VASP).15−18 The Perdew−Burke−Ernzerhof (PBE) functional19 was used in the DFT+U approach of Dudarev et al.20 with a U value of 5 eV. We do not expect that the results presented here depend critically on the choice of the U value. This is due to the fact that the U value mainly affects energy differences in situations where the number of occupied f electrons changes, as for instance, in the creation of an oxygen vacancy. In the expression for the formation energies calculated in this work, the number of f electrons is conserved and consequently we expect a very weak U dependence. In a recent publication,7 we confirmed that this is indeed the case for the fluorine dimer. The fluorine dimer formation energy was calculated to be −0.09 with U = 3 eV and −0.10 with U = 5 eV, respectively. A cutoff of 400 eV was imposed for the plane-wave basis set describing the valence electrons (12 electrons for cerium, 6 for oxygen, and 1 for hydrogen), while the core electrons were treated with the projector augmented wave method (PAW) method developed by Blöchl.21 Aspherical contributions from the gradient corrections inside the PAW spheres were included. The (111) surface was modeled with full 3D periodicity and a lamellar construction with slabs that were three O−Ce−O triple layers thick, with ∼15 Å thick vacuum gaps between them. A p(6 × 6) surface supercell was used, and the cell parameters were kept fixed at the PBE+U optimized bulk value of 5.50 Å. The F defects were introduced at one side of the slab, and the positions of all atoms in the slabs were optimized until the forces on them were less than 0.01 eV/Å. In all cases, the Brillouin zone was sampled at the Γ-point. Systems involving localized electrons in the form of Ce3+ ions need to be treated with extra care; here possible problems with multiple local electronic minima were solved using the so-called occupation matrix control, originally designed for UO2,22 but later also for ceria.23 We used the implementation from ref 22.

ϕ(a , b) = Ca,b for nearest neighbors =−

qaqb rab

otherwise (2)

Here qa is the charge of species a and rab is the distance between a and b; charges are such that qF = −qCe3+ = q, with q > 0. The four parameters {CF,F, CCe3+,Ce3+, CF,Ce3+, q} in this model were adjusted to reproduce a set of 78 cluster formation energies obtained from DFT calculations. Full Ewald summations were used during the fitting process as well as in the production runs. The physical basis for the form of the interaction model is rather straightforward. Stoichiometric ceria is made up of Ce4+ and O2− ions. Substitution of an O by an F produces an F− ion and a Ce3+ ion, that is, a charge change of +1 and −1, respectively, with respect to the ceria lattice. At large distances between the various substituents we expect the main contribution to interaction energy differences to be pairwise additive and Coulombic in form. However, the F− and Ce3+ substituents also induce lattice and electronic distortions in their respective immediate surroundings. As long as the defect ions remain sufficiently far away from each other, the energy contributions from such relaxations can be expected to be appreciable but equal to those around isolated defects and thus constant and not contributing to the energy change when the defects move around. If we bring, say, two F− substituents at the ceria surface into nearest-neighbor position to each other, however, then the local lattice and electron relaxations will be very different from a superposition of the two individual relaxations that we had before. As a result, besides the Coulomb contribution, we now have to take into account the additional energy contribution due to the nonadditivity of these relaxations. Because even the concept of Coulomb interaction between individual charges may have lost its conceptual meaning at this distance, we chose not to partition the energy into two contributions but just use one for each pair of species, viz. F−−F−, Ce3+−Ce3+, and F−−Ce3+. This discussion of our rationale behind the form of our model also gives some insight into the limitations of the model. We have explored several other interaction models of varying complexity, also including additional fitted parameters for larger separations. However, we found that most performed worse at reproducing the DFT clustering energies and none performed significantly better. Given the limitations of a pairwise interaction model and the simplicity of the specific form we have chosen, the interaction model reproduces the quantummechanical energies surprisingly well, as will be shown in the Results section. DFT Calculations. We created a training set of defect cluster formation energies from quantum-mechanical calculations of the DFT+U type for subsequent use in the fitting of a simple classical interaction model (as described above). Our definition of the quantum-mechanical formation energy of a defect cluster (F−)n(Ce3+)n is EDFT = Etot(slab with (F−)n (Ce3 +)n defect cluster) − nEtot(slab with one isolated F−Ce3 + defect) + (n − 1)Etot(stoichiometric slab)

(3)

Figure 4. Average fraction of F− monomers relative to that of a random distribution, as a function of temperature, for both the “standard” Monte Carlo simulation (solid line) and the “minimumenergy cerium configuration” Monte Carlo simulation described in the Supporting Information (filled circles).

The DFT+U formation energies were calculated for nearestneighbor (NN) F dimer defect clusters and for all symmetrically inequivalent connected F trimer clusters. The clusters are shown in Figure 1 together with the associated Ce3+ ions. The D


Article


Figure 5. Results from Monte Carlo simulations at 600 K and with 15% F concentration. A large number of snapshots from the simulations were analysed with respect to the number, NF, of F monomers present, and their relative probabilities were calculated. This resulted in an approximately Gaussian probability distribution function with average number of F monomers, ⟨NF⟩, equal to 85 and standard deviation of σ = 7.5. The eight snapshots shown were chosen from those for which NF ∈ (NF − σ, NF + σ).

■

RESULTS DFT Results and Interaction Model. Analysis of the resulting DFT energies shows that, for all F cluster configurations, the most stable Ce3+ configuration is that which maximizes the F−−Ce3+ coordination. Thus the Ce3+ ions act as the “glue” that binds the F− together, and the cluster formation energy depends rather significantly on the positions of the Ce3+ ions. For example, the DFT-calculated formation energies actually become positive (unstable) if the Ce3+ ions are located around the edges of the clusters, as might be expected given the repulsive interaction between the F− ions themselves. The four parameters in our potential energy model were determined by minimizing the expression

∑ [EiDFT − Eimodel]2

Figure 6. Results from the DFT calculations and Monte Carlo simulations of ceria(111) with surface OH groups. (a) Correlation plot for the formation energies per OH of the 12 defect clusters in the training set. The DFT energies from the training set are on the x axis and the energies of the fitted model function (from expressions 1, 2, and 5 in the text) are on the y axis. For the reader’s convenience, two energy scales are given in the Figure. (b) Comparison of a snapshot from the Monte Carlo simulation of OH at a temperature of 600 K and a defect concentration of 15%, with the snapshot from the corresponding F simulation, copied from Figure 3.

(4)

i

with the sum running through our set of 78 defect configurations. EDFT is the quantum-mechanical defect cluster i formation energy of cluster i as defined in equation 3, and Eimodel is the corresponding energy calculated from the interaction model. A complete overview of the DFT clustering energies is given in the Supporting Information. The final parameters resulting from our fitting procedure are C F,F = 0.299 eV C F,Ce3+ = −0.370 eV

clusters. The agreement between model and reference DFT energies in Figure 2 is seen to be good. This gives us confidence to proceed to the Monte Carlo simulations. Clustering Behavior. Figure 3 presents results from Monte Carlo simulations carried out at different temperatures (300, 600, and 900 K) and different fluorine defect concentrations (5, 10, and 15% out of the total number of anionic sites in the top layer). For comparison, similar results for random distributions with the same defect concentrations are shown in the top row of the Figure. For each simulation, the left panel gives the fraction of all F− ions that take part in clusters ranging from monomers to hexamers. A fluorine cluster is defined as a group of F− ions in which any F pair may be connected by a path of nearest neighbor steps within the cluster. Furthermore, all clusters of sizes larger than dimers are classified as being either open, when no three among the F− ions in the cluster form a triangle, or

CCe3+,Ce3+ = 0.220 eV q = 0.276 e

(5)

which yields Φmodel in electronvolts. We note that the fitting procedure leads to a noninteger value for the charge. Alternatively, the charges may be set to ±1, and one can use a relative (static) dielectric constant instead. This would give a value for the latter of 1/(0.276)2 = 13.2, which is not unreasonable given that the bulk dielectric constant of ceria is on the order of 20 and that an effective dielectric constant at the surface should be lower due to the smaller quantity of material surrounding the charged species. The quality of our fitted model can be appreciated from Figure 2, where we have plotted the defect formation energy per flourine with our model (Eimodel/n) against the corresponding quantum-mechanical energy (EDFT /n) of our set of 78 fluorine dimer and trimer i E


Article

The Journal of Physical Chemistry C compact otherwise. In each bin, the contributions of open and compact clusters are indicated by white and dashed markings, respectively, as indicated for the (600 K, 5% F) case. The representative snapshot to the right of each histogram, was chosen to be that one among all the Monte Carlo frames whose cluster distribution was the most similar to the corresponding run time average shown to the left of it, that is, such that the sum of squared deviations between the two is minimal. The cluster histogram for the 10% concentration case at 300 K is very flat, meaning that approximately equally many fluorines belong to clusters of all sizes up to a maximum of ∼25 (not shown). Most of the larger clusters are rather compact or consist of two or more connected compact pieces. At 600 K, the distribution has shifted toward smaller clusters, and the trimers mainly adopt open, linear structures instead of the compact triangular structures. The tetramers are also more open/linear than at the lower temperature. At 900 K, the same system is even richer in small clusters and isolated fluorines, with hardly any compact clusters remaining. The same development with temperature occurs for the two other F concentrations. The top row of Figure 3 shows histograms and representative snapshots of randomly produced F distributions (as usual subject to a constraint of 0 or 1 occupancy of the sites). Clearly, for all coverages, with increasing temperatures the histograms and snapshots for our model systems begin to resemble those of the random distributions. However, for all F concentrations in Figure 3, the number of isolated defects is larger at 900 K than for the random distribution. This result might appear counterintuitive, and to explore this behavior further we performed additional Monte Carlo simulations for the 15% case at higher temperatures, up to almost 10 000 K. The result is given in Figure 4. The solid curve in this Figure shows the average number of isolated fluorines (monomers) at the ceria surface as a function of temperature, plotted relative to the number of monomers in the random distribution (dashed line). While this ratio is 1 beyond Ttrans = 700 K and then finally decays slowly toward one, as the temperature increases further. The presence of more single defects than in a random distribution clearly indicates that the effective interaction between them has become repulsive. An intuitive explanation of this behavior may be obtained by invoking the “glue” concept introduced in the DFT Results and Interaction Model sub-section. The system consists of two sublattices, each populated with mutually repulsive particles, while the particles on the two different sublattices interact attractively. This latter attraction produces the glue that creates the clusters on the individual lattices. With increasing temperature the thermal agitation becomes strong enough to overcome the attractions between the sublattices and make the Ce3+ leave the energetically most stable positions nearest to the fluorines, exposing the F ions to each other. From a thermodynamic point of view, the attractive−repulsive transition is induced by the increasing importance of the entropy of the Ce3+. This picture is corroborated by the fact that simulations with the Ce3+ restricted to their energetically most stable positions, do not have this transition. Those simulations are illustrated by the red dots in Figure 4, and described in greater detail in the Supporting Information. Comparison with Experiment. We next analyze the case of 15% coverage at a temperature of 600 K in some more detail.

Interestingly, even though compact clusters are energetically more stable than open clusters, we find many linear or hooked structures under these conditions, sometimes even consisting of as many as five to six F atoms. Incidentally, note that this happens despite the fact that our interaction model overestimates the binding energy of the triangular configurations (see the leftmost point in Figure 2). The abundance of linear structures may not be very clear from the snapshots in Figure 3 at first sight, but we stress that this to a large extent is due to the selection protocol used in that Figure. Recall that snapshots were chosen such that their cluster size distributions were as close as possible to the run-time-averaged histogram. This tends to hide “extreme” features that otherwise would have immediately caught the human eye in a visual inspection. We therefore present in Figure 5 a few additional snapshots from the Monte Carlo simulation at 600 K and with a 15% concentration to show more clearly the abundance of linear and hooked structures. They were randomly chosen from a subset of all snapshots whose fraction of F monomers was within one standard deviation of the run-time average. These snapshots clearly show that even when being selected to have many F monomers present, still many extended and hooked linear structures occur, similar to the ones shown in the paper by Esch et al. It is important to note that our linear structures start to appear at temperatures between 600 and 900 K, within the range of temperatures probed by the Esch experiment.4 Finally, we note that OH groups give very similar STM signatures to oxygen vacancies and F impurities9 and therefore are alternative candidates for the features seen in STM pictures. To explore this possibility, we have performed preliminary simulations with surface hydroxyls following a similar approach to that described above for F impurities. The specific parameters in our potential energy model for OH groups were fitted to a data set of 12 OH clusters (see Supporting Information for more details). For consistency, the charge, q, and the CCe3+,Ce3+ parameters were kept fixed at the optimized values from the fluorine model described above. The values for the fitted OH parameters were COH,OH = 0.295 and COH,Ce3+ = −0.314 eV. Figure 6a shows a correlation plot demonstrating the quality of the fitted parameters, while Figure 6b compares a snapshot from the OH simulations at a defect concentration of 15% and a temperature of 600 K with a corresponding snapshot from the fluorine simulations. Clearly, OH groups form far fewer and smaller nearest neighbor clusters than F, or than the random distribution for that matter. Instead, OH groups tend to form configurations where they are located as next-nearest neighbors. Consequently, our hydroxyl simulations predict that surface hydroxyls are unlikely to be the main impurities seen in the experiments of Esch et al.

■

SUMMARY AND CONCLUDING REMARKS We have studied the distribution of fluorine substitutional impurities in the top layer of a CeO2(111) surface, using a combination of DFT+U and lattice Monte Carlo simulations, which take into account the effects of configurational entropy of both the fluorines themselves and the Ce3+ ions. Using DFT, we calculated the formation energies of 78 surface clusters of fluorine defects, up to trimers in size, with the associated Ce3+ ions within nearest-neighbor distance of at least one F− ion. We found that in the most stable configurations the Ce3+ maximize their coordination to the F−, acting in some sense like glue holding them together. The energy differences between different cluster shapes and F


Article

The Journal of Physical Chemistry C between different Ce3+ configurations are all rather small, namely, on the order of a few 100 to 1000 K, which implies that entropic effects should play a major role in determining the F distribution at all but very low temperatures. To be able to perform Monte Carlo simulations, we constructed a pairwise additive force-field model for systems of fluorines and Ce3+ ions, taking into account the electrostatic interactions between the various species as well as interactions between their strain fields at short distances. The force-field model parameters were adjusted to obtain the best agreement between model energies and DFT energies for the set of 78 configurations. Given its simplicity, the interaction model performs surprisingly well in reproducing the DFT data. An interesting result from the Monte Carlo simulations relates to the number of isolated fluorines, which actually overshoots the number expected from a random distribution at high temperatures. By carrying out alternative Monte Carlo simulations (Supporting Information), in which the Ce3+ always seek their most energetically favorable locations regardless of the fluorine temperature, we found that the larger-than-random number of “isolated” monomers is a direct result of the Ce3+ configurational entropy. Overall, the Monte Carlo simulations showed that, for all defect concentrations considered, the distribution of fluorines shifts toward smaller and more open clusters as the temperature increases. In the temperature range 600−900 K the number of open (linear or bent) trimer or larger clusters found is still appreciable, namely, ∼20% at a defect concentration of 15%. Some eye-catching sample structures were shown in Figure 5. We conclude from our analysis that the abundance of open clusters observed in experiments is mainly a result of entropy rather than energetics.

■

Julian Gale (Curtin University), who readily modified the GULP code for our purposes. J.K. also acknowledges Dr. Meysam Pazoki for useful comments.

■

(1) Nörenberg, H.; Briggs, G. A. D. Defect Structure of Nonstoichiometric CeO2(111) Surfaces Studied by Scanning Tunneling Microscopy. Phys. Rev. Lett. 1997, 79, 4222−4225. (2) Nörenberg, H.; Briggs, G. Defect Formation on CeO2 Surfaces After Annealing Studied by STM. Surf. Sci. 1999, 424, L352−L355. (3) Namai, Y.; Fukui, K.-i.; Iwasawa, Y. Atom-Resolved Noncontact Atomic Force Microscopic Observations of CeO2(111) Surfaces with Different Oxidation States: Surface Structure and Behavior of Surface Oxygen Atoms. J. Phys. Chem. B 2003, 107, 11666−11673. (4) Esch, F.; Fabris, S.; Zhou, L.; Montini, T.; Africh, C.; Fornasiero, P.; Comelli, G.; Rosei, R. Electron Localization Determines Defect Formation on Ceria Substrates. Science 2005, 309, 752. (5) Conesa, J. Computer Modeling of Surfaces and Defects on Cerium Dioxide. Surf. Sci. 1995, 339, 337−352. (6) Zhang, C.; Michaelides, A.; King, D. A.; Jenkins, S. J. Oxygen vacancy clusters on ceria: Decisive role of cerium f electrons. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 79, 075433. (7) Wolf, M. J.; Kullgren, J.; Broqvist, P.; Hermansson, K. Fluorine impurities at CeO2(111): Effects on oxygen vacancy formation, molecular adsorption, and surface re-oxidation. J. Chem. Phys. 2017, 146, 044703. (8) Conesa, J. C. Surface Anion Vacancies on Ceria: Quantum Modelling of Mutual Interactions and Oxygen Adsorption. Catal. Today 2009, 143, 315−325. (9) Kullgren, J.; Wolf, M. J.; Castleton, C. W. M.; Mitev, P.; Briels, W. J.; Hermansson, K. Oxygen Vacancies versus Fluorine at CeO2(111): A Case of Mistaken Identity? Phys. Rev. Lett. 2014, 112, 156102. (10) Wu, X.-P.; Gong, X.-Q. Clustering of Oxygen Vacancies at CeO2: Critical Role of Hydroxyls. Phys. Rev. Lett. 2016, 116, 086102. (11) Wolf, M. J.; Kullgren, J.; Hermansson, K. Comment on “Clustering of Oxygen Vacancies at CeO2: Critical Role of Hydroxyls. Phys. Rev. Lett. 2016, 117, 279601. (12) Grieshammer, S.; Grope, B. O. H.; Koettgen, J.; Martin, M. A combined DFT + U and Monte Carlo study on rare earth doped ceria. Phys. Chem. Chem. Phys. 2014, 16, 9974−9986. (13) Purton, J. A.; Allan, N. L.; Gunn, D. S. D. Simulations of Doped CeO2 at Finite Dopant Concentrations. Solid State Ionics 2017, 299, 32−37. (14) Gale, J. D.; Rohl, A. L. The General Utility Lattice Program (GULP). Mol. Simul. 2003, 29, 291−341. (15) Kresse, G.; Hafner, J. Ab Initio Molecular-Dynamics Simulation of the Liquid-Metal−Amorphous-Semiconductor Transition in Germanium. Phys. Rev. B: Condens. Matter Mater. Phys. 1994, 49, 14251−14269. (16) Kresse, G.; Hafner, J. Ab Initio Molecular Dynamics for Liquid Metals. Phys. Rev. B: Condens. Matter Mater. Phys. 1993, 47, 558−561. (17) Kresse, G.; Furthmüller, J. Efficient Iterative Schemes for ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 11169−11186. (18) Kresse, G.; Furthmüller, J. Efficiency of Ab-Initio Total Energy Calculations for Metals and Semiconductors Using a Plane-Wave Basis Set. Comput. Mater. Sci. 1996, 6, 15−50. (19) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (20) 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: Condens. Matter Mater. Phys. 1998, 57, 1505−1509. (21) Blöchl, P. E. Projector Augmented-Wave Method. Phys. Rev. B: Condens. Matter Mater. Phys. 1994, 50, 17953−17979. (22) Dorado, B.; Amadon, B.; Freyss, M.; Bertolus, M. DFT+U Calculations of the Ground State and Metastable States of Uranium Dioxide. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 79, 235125.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.7b00299. Table S1. Binding energies for the surface fluorine (OH) dimer. Table S2. Binding energies for the surface fluorine (OH) linear trimer. Figure S1. Schematic plan view of the uppermost triple layer of the CeO2(111) surface facet. Table S3. Binding energies for the surface fluorine hooked trimer. Table S4. Binding energies for the surface fluorine (OH)“triangle 1” trimer. Table S5. Binding energies for the surface fluorine (OH)“triangle 2” trimer. (PDF)

■

REFERENCES

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Jolla Kullgren: 0000-0003-3570-0050 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported the Swedish Research Council (VR). Some of the simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at UPPMAX and NSC. We acknowledge the National Strategic e-Science program eSSENCE. Special thanks go to Professor G


Article

The Journal of Physical Chemistry C (23) Allen, J. P.; Watson, G. W. Occupation Matrix Control of d- and f-electron Localisations Using DFT+U. Phys. Chem. Chem. Phys. 2014, 16, 21016−21031.

H


DFT-based Monte Carlo Simulations of Impurity ... - ACS Publications

Recommend Documents