This item is only available as the following downloads:
1 REACTIVE VARIABLE CHARGE POTENTIAL DEVELOPMENT AND ATOMISTIC S IMULATIONS OF SURFACE AND HETER O GENEOUS INTERFAC IAL INTERACTIONS By YU TING CHENG A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2013
2 2013 Yu Ting Cheng
3 To my family with love
4 ACKNOWLEDGMENTS First, I would like to show my gratitude to my advisor Prof. Susan B. Sinnott for h er kind instruction and all efforts to guide me during my Ph.D. study Her continuous support s suggestions and helpful discussion s ha ve helped me pass through all difficulties during my research experience s I would also like to thank Prof. Simon R. Phillpot for his invaluable advices on the studies in this dissertation and to extend my sincere appreciation to my committee members Prof. Paul H. Holloway, Prof. Scott S. Perry and Prof. Sergey Vasenkov for their supports and comments Great thanks to all former and present members in the Computational Materials Science Focus Group for their fr iendship and continu ous support s. I also want to specially thank Dr s Tao Liang, Tzu Ray Shan and Bryce Devine for their mentor and help on my researches, Te Yu Kao, Ya Chiao Chang, P a ng Wei Liu, Hsiu Wen Tsai, Aiden Liu, Mong Jen C hen, Hsin Yin Tsai, Yi Chung Wang, Ya Wen Yeh, Chung Hsin Fan, Chao Hsuan Shen, Kun Ta Tsai, Tsu Wu Chiang and Tzung Han Lai for their constant support and encourage inside and outside the school. Most importantly my beloved family give s me unconditional love, patien ce and support Without them, none of this would have been possible Finally I want to express my gratitude to my lovely wife, Chih Yun (Angela) Huang, for h er never ending support The dissertation is dedicated to her.
5 TABLE OF CONTENTS page ACKNOWLEDGMENTS ................................ ................................ ................................ .. 4 LIST OF TABLES ................................ ................................ ................................ ............ 8 LIST OF FIGURES ................................ ................................ ................................ ........ 10 ABSTRACT ................................ ................................ ................................ ................... 12 CHAPTER 1 INTRODUCTION ................................ ................................ ................................ .... 14 1.1 General Introduction ................................ ................................ ......................... 14 1.2 Integration of Computational and Experimental Methods ................................ 15 1.3 Outline ................................ ................................ ................................ .............. 18 2 COMPUTATIONAL METHODS ................................ ................................ .............. 25 2.1 Computational Methods ................................ ................................ .................... 25 2.2 Density Functional Theory ................................ ................................ ................ 25 2.3 Classical Empirical Methods ................................ ................................ ............. 27 2.3.1 Interatomic Potentials ................................ ................................ .............. 27 2.3.2 Molecular Dynamics (MD) Simulations ................................ .................... 28 2.3.3 Ensembles ................................ ................................ ............................... 29 2.3.4 Periodic Boundary Conditions and Cutoff Radius ................................ .... 30 2.4 Variable Charge Reactive Potentials: COMB Potentials ................................ ... 30 2.4.1 General Formalism of the Third Generatio n COMB Potentials ................ 30 2.4.2 Electrostatic Energies ................................ ................................ .............. 31 2.4.3 Short Range interactions ................................ ................................ ......... 34 2.4.4 v an d er Waals Interaction s ................................ ................................ ...... 36 2.4.5 Correction Terms ................................ ................................ ..................... 37 2.4.6 Dynamic Charge Equilibration Method ................................ .................... 37 2.5 Parameterization of COMB Potentials ................................ .............................. 39 2.5.1 Creating the Database ................................ ................................ ............. 40 2.5.2 Parameterization of Pure Systems ................................ .......................... 41 2.5.3 Parameterization of Binary Systems ................................ ........................ 41 2. 6 Transition State Search: Dimer Method ................................ ............................ 42 2. 7 Summary ................................ ................................ ................................ .......... 43 3 DEVELOPMENT OF COMB3 CU/ZNO POTENTIAL AND ITS APPLICATION ON CU CLUSTERS DEPOSITION ON ZNO SURFACE ................................ ........ 50 3.1 Potential development and Properties of Cu/Z nO ................................ ............. 52 3.1.1 Fitting Procedure ................................ ................................ ..................... 52 3.1.2 Predicted Properties of ZnO polymorphs ................................ ................. 53
6 3.1.3 Predicted Properties of CuZn alloys ................................ ........................ 54 3.2 Interaction of Cu on ZnO ................................ ................................ ....... 54 3.2.1 Adsorption Energies of Cu Adatom on ZnO ................................ 54 3.2.2 Migration barriers of Cu Adatom on ZnO ................................ ..... 55 3.3 Cu Clusters Deposition on ZnO ................................ ............................. 56 3.3.1 Details of the MD Simulation of Cluster Deposition and Relaxation ........ 56 3.3.2 Surface Morphology and Growth Mode ................................ ................... 57 3.3.3 Effect of Incident Energy ................................ ................................ ......... 59 3.3.4 Effect of Surface Temperature ................................ ................................ 60 3.4 Comparison of the Growth of Cu on Cu(111) and ZnO ......................... 61 3.5 Summary ................................ ................................ ................................ .......... 62 4 DEVELOPMENT OF COMB3 TI/TIO 2 POTENTIAL AND ITS APPLICATONS ON CU CLUSTERS ON TIO 2 (110) SURFACES ................................ ..................... 76 4.1 Potential Development and Properties of Ti and TiO 2 ................................ ....... 77 4.1.1 Fitting Procedure ................................ ................................ ..................... 77 4.1.2 Predicted Ti Properties ................................ ................................ ............ 78 4.1.3 Properties of Rutile TiO 2 and Phase Stability ................................ .......... 80 4.1.4. TiO 2 Low Miller Index Surfaces ................................ .............................. 81 4.2 Heterogeneous Interfacial System: Cu n /TiO 2 (110) ................................ ............ 83 4.2.1 Cu Clusters on Stoichiometric TiO 2 (110) ................................ ................. 83 4.2.2 Cu 147 Clusters on Non Stoichiometric TiO 2 (110) ................................ ..... 84 4.3 Summary ................................ ................................ ................................ .......... 85 5 DEVELOPMENT OF COMB3 TI/TIN POTENTIAL AND ITS APPLICATONS ON TI/TIN INTERAFACES AND THE DEPOSTION OF ATOMIC OXYGEN ON TIN(001) SURFACE ................................ ................................ ................................ 95 5.1 Potential development and Properties of Ti, N 2 and TiN ................................ .. 96 5.1.1 Fitting Procedure ................................ ................................ ..................... 96 5.1.2 Fitted Properties of Ti metal and N 2 molecule ................................ ......... 97 5.1.3 Fitted and Predicted TiN Phase Stability and Properties ......................... 98 5.1.4 Predicted Defect Formation Energies and Surface Energies of rocksalt TiN ................................ ................................ ................................ ... 99 5.2 Predicted Dynamical Results ................................ ................................ .......... 101 5.3 Applications of COMB 3 TiN Potential ................................ ............................. 102 5.3.1 Adhesion of fcc Ti/TiN(001) Interfaces ................................ .................. 102 5.3.2 Atomic Oxygen deposition on TiN(001) surface ................................ .... 104 5.4 Summary ................................ ................................ ................................ ........ 105 6 FIRST PRINCIPLES STUDY OF EARLY STAGE OXIDATION OF FACETED AlN SURFACES ................................ ................................ ................................ ... 113 6.1 Computational Details ................................ ................................ ..................... 113 6.2 Adsorption of Atomic and Molecular Oxygen on the A lN surface .................... 115
7 6.3 Summary ................................ ................................ ................................ ........ 116 7 CONCLUSIONS ................................ ................................ ................................ ... 120 APPENDIX : POTEN TIAL PARAMETERS FOR COMB POTENTIALS DEVELOPED IN THIS WORK ................................ ................................ ................................ ..... 123 REFERENCE LIST ................................ ................................ ................................ ...... 128 BIOGRAPHICAL SKETCH ................................ ................................ .......................... 141
8 LIST OF TABLES Table page 1 1 Characteristics of physical adsorption and chemisorption. ................................ 20 2 1 The com mon used interatomic potentials ................................ ........................... 45 2 2 List of fitting parameters for the pure systems. ................................ ................... 46 2 3 List of fitting parameters for the charge dependent energy terms. ..................... 46 2 4 List of fitting parameters for the binary systems. ................................ ................ 46 3 1 Comparison of the properties of ZnO wurtzite phase predicted by the COMB2 and COMB3 potentials with values obtained from experiments and DFT calculations. ................................ ................................ ................................ ........ 64 3 2 Properties of Cu Zn alloys predicted by COMB3 compared to those from e xperiments and DFT calculations ................................ ................................ ..... 64 3 3 The adsorption energy (eV) of a single Cu adatom at the three different adsorption sites on the ZnO surface as predicted by COMB3 and DFT calculations ................................ ................................ ................................ ......... 65 3 4 Predicted migration barriers of a Cu adatom on the ZnO surface found using the dimer method with COMB 3 and compared with DFT results .... 65 4 1 Properties of Ti reproduced or predicted by the COMB3 potential developed in this work, and compared to experimental results, DFT calculations and the MEAM potential. ................................ ................................ ................................ 87 4 2 Properties of the rutile, anatase, and brookite phases calculated by the COMB potential. The properties are compared to experiments, DFT calculations and three other empirical potentials. ................................ ............... 88 4 3 Cohesive energ ies of rutile titania and that relative to rutile titania for other polymorphs calculated from COMB potential in comparison with DFT calculations and other empirical potentials ................................ ........................ 89 4 4 Low index surface energies of rutile TiO 2 from the ab initio approaches and empirical potentials. ................................ ................................ ............................ 90 4 5 Charge variation (Q Q bulk ) of surface atoms of TiO 2 surfaces ............................. 91 5 1 COMB predictions compared to published experimental results, DFT calculations and MEAM potentials for the lattice constants a 0 and c 0
9 enthalpy of formation, cohesive energy, bulk modulus, shear modulus and el astic constants of the rocksalt TiN and rutile Ti 2 N ................................ .......... 107 5 2 Comparison o f formation energies of poin t defects in NaCl typed TiN crystal calculated from first principles calculations and COMB potentials. ................... 108 5 3 Relaxed surface energy of TiN (100), (110) and (111) surface predicted by COMB potentials, in comparison with the predictions from published first princip les calculation and MEAM potentials. ................................ .................... 108 6 1 Comparison of the computed bulk properties to the experimental data. ........... 117 6 2 Predicted adsorption energies (E ads ) of atomic oxygen on the flat and stepped AlN surfaces. ................................ ................................ ........... 117 A 1 Atomic and electrostatic parameters of pure elements (O, Cu, N, Ti, and Zn) for COMB3 potentials. ................................ ................................ ...................... 124 A 2 Bond typed parameters of pure element (O, Cu, N, Ti, and Zn) for COMB3 potentials. ................................ ................................ ................................ ......... 125 A 3 Bond typed parameters of Cu/ZnO system for COMB potential developed in this work. ................................ ................................ ................................ .......... 126 A 4 Bond typed parameters of TiN and TiO 2 systems for COMB3 potentials developed in this work. ................................ ................................ ..................... 127
10 LIST OF FIGURES Figure page 1 1 Initial oxidation processes of bare metal surfaces. ................................ ............. 21 1 2 Surface processes of metal thin films or clusters growth on metal oxide substrate ................................ ................................ ................................ ............. 22 1 3 Hierarchy of multiscale modeling ................................ ................................ ........ 23 1 4 Snapshot of a molecular dynamics simulations of water, carbon monoxide, and carbon dioxide molecules interacting with a copper catalyst supp orted by a zinc oxide substrate ................................ ................................ ......................... 24 2 1 Lennard Jones (LJ) potential ................................ ................................ ............. 47 2 2 Schematic representation of periodic bounda ry conditions in MD simulations ... 48 2 3 The two migration mechanisms of single Cu atom on Cu(100). ......................... 49 3 1 Top view of three different adsorption sites of Cu atoms on the ZnO surface ................................ ................................ ................................ ................ 66 3 2 The binding energy (eV) of Cu n (n = 2 to 5) clusters on ZnO predicted by DFT calculations and COMB potential. ................................ .......................... 66 3 3 Relaxed structure of Cu 5 clusters supported by the ZnO surface. .......... 67 3 4 Simulation model for deposition of Cu clusters on the ZnO surface. ....... 68 3 5 MD simulation of the morphology of Cu clusters following the deposition of 1ML of Cu atoms. ................................ ................................ ............................... 69 3 6 Snapshots of just the Cu atoms at the indicated coverage (0.2 1ML) following deposition on ZnO ................................ ................................ .... 69 3 7 The deposition processes of Cu clusters landing on ZnO surface ...................... 70 3 8 Schematic of Ehrlich Schwoebel barrier ( ES ) ................................ ................. 70 3 9 Top down snapshots from MD simulation of the surface morphology of 1ML of Cu deposited on ZnO with different incident energies. ........................ 71 3 10 Top down snapshots from MD simulation of the morphology of 1ML Cu deposited on ZnO maintained at the indicated temperatures .................. 7 2
11 3 11 Snapotshot from an MD simulation of the morphology of Cu clusters on ZnO at the indicated coverage. ................................ ................................ ...... 73 3 12 Snapotshot from an MD simulation of the morphology of Cu clusters on Cu(111) at the indicated coverage ................................ ................................ ...... 74 3 13 Top down snapshots from the MD simulations of 2ML Cu deposited on Cu(111) and ZnO ................................ ................................ .................... 75 4 1 Model of the rutile TiO 2 surfaces ................................ ................................ ....... 92 4 2 Optimized geometries of Cu 13 Cu 38 Cu 55 and Cu 147 clusters. .......................... 93 4 3 Charge distribution of relaxed Cu clusters on TiO 2 (110) ................................ .... 93 4 4 Charge distribution of relaxed Cu 55 clusters on the stoichiometric and defected TiO 2 (110) surfaces ................................ ................................ ............... 94 5 1 The unit cell of fitting structures for TiN and Ti 2 N systems in this study ........... 109 5 2 Schematic of layers sequence of the TiN(001), (110) (111) surfaces .............. 109 5 3 The crystal structures of rocksalt T iN as a function of temperature .................. 110 5 4 The interfacial structures of Ti(001)/TiN(001) ................................ ................... 111 5 5 Snapshots of oxygen atoms de posited on the TiN(001) surface ...................... 112 5 6 Snapshots from the deposition of single oxygen molecule on the TiN(001) surface ................................ ................................ ................................ .............. 112 6 1 Stable adsorption sites on the flat and st epped AlN surfaces ............... 118 6 2 The energy chage of the adsorption and dissociation of O 2 on the flat and stepped AlN surface ................................ ................................ .............. 119
12 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy REACTIVE VAR IABLE CHARGE POTENTIAL DEVELOPMENT AND ATOMISTIC SIMULATION S OF SURFACE AND HETEROGENEOUS INTERFACIAL INTERACTIONS By Yu Ting Cheng December 2013 Chair: Susan B. Sinnott Major: Materials Science and Engineering A dvancement s in computational resources ha ve enabled the application of computational methods in a manner that complement s experimental measurements and ultimately captur e the process structure property relationship s of material systems Among various computational methodologies the classical empirical methods (also known as interatomic potentials) has been developed and used in molecular dynamics (MD) simulations Such simulations are able to model systems at atomic and nanometer scale and dynamical process es that include the effects of temperatures and pressure which is beyond the reach of quantum based approaches. Because the interatomic potential is the main gradient in MD simulation s for calculating energies and forces, it is therefore important that the interatomic potential is accurate and transferable in order to correctly describe the chemical and physics properties of materials under a variety of conditions. In addition, to ensure the wide st possible application, the interatomic potential should be ca pable of simulating system s with different types of chemical bonding present
13 Recently, the charge optimized many body (COMB) potentials were developed to tackle this challenge and parameterized to include a wide range of elements and compounds, including heterogeneous systems. I n this dissertation, the main objective is to apply and parameterize the third generation COMB potential for several specific heterogeneous systems and use them to investigate the underlying physics of interfac ial and surface phenomena that are critical for the success of many applications P otentials were developed for metallic Zn and Ti, and for CuZn, CuTi, ZnO, TiO 2 and TiN Through MD simulations COMB potentials have been successfully applied to characterize the growth mode s of Cu on Cu (111) and ZnO and clearly elucidate the ways in which Cu grow th transitions from layer by layer to three dimensional as Cu coverage increases. For the Cu/TiO 2 system, a n enhanced bonding between Cu clusters and the oxidized TiO 2 (110) surface was predicted T he adsorption mechanism was predicted to be th e formation of a metal ox ygen bond between Cu and O. In addition, the characterization of the adsorption of O and O 2 on the TiN(001) surface showed that the COMB potentials for Ti/TiN/TiO 2 and N/NO systems were well suited to the oxidation study. As a result of the work carried out and described in this dissertation, COMB potentials can be considered to be potent tools from an engineering perspective to provide useful guidance for the interpretation of experimental data and from a fundamental scien ce perspective to facilitate materials development.
14 CHAPTER 1 INTRODUCTION 1.1 General Introduction Processes at surface are critical for the success of many applications 1 5 For example surface oxidation in particular at high temperature s is inevitable for many systems w here the application takes place under ambient conditions. The process of surface oxidation may result in destructive outcomes such as the degradation of thermal barrier coatings 6 7 or beneficial outcomes such as heterogeneous catalysts. 8 9 With recent advances in material characterization techniques, such as in situ transmission electron microscopy (TEM) Auger electron spectroscopy (AES), and X ray photoelectron spectroscopy (XPS), experimental methods are able to provide atomic resolution imag es of surfaces and acc ompanying information regarding the details of their structure, electronic structure, and composition. As a result, s everal mechanisms of the initial oxidation process on bare metal surfaces have been determined, 10 12 as shown in Figure 1 1 In principle when oxygen molecules impinge on material surfaces, two adsorption process es may take place depending on type of interaction that are prevalent. 13 One is physi sorption where the interaction s are dominated by van der Waals forces. Because of the weak interaction the molecule has a high mobility on surfaces and th is adsorption is reversible which leads to high de sorption rates In contrast chemisorption involves the transfer and sharing of electrons which leads to c hemical bond formation Therefore chemisorption results in stronger molecule surface interactions th a n physi sorption. In addition, chemisorption is sensitive to details of the s urface structure and its symmetry. The detaile d characteristics of physi sorption and
15 chemisorption are listed in Table 1 1. Subsequently, the next process is metal oxide formation, for which two possible mechanism s have been proposed. First, the dissociative oxygen atoms penetrate into the subsurface and form an oxide within the sub strate Second, the substrate surface atoms exchange with adsorped oxygen atoms and form an oxide film at the surface. Many of the a bove describe d processes at sur face are also important in heterogeneous catal ysis using metal/metal oxide systems Reactions where such catalysts are important include the synthesis of methanol via hydrogenation of CO and CO 2 14 water gas shift reaction for removal of CO 15 16 and the production of hydrogen by steam reforming of m ethanol 17 Despite the technological importance of these r eactions, in many instances the nature of the active sites and their relative stability is still not well established. Factors such as surface geometry, the interface between the metal and oxide, and the size and distribution of metal clusters strongly con trol the reactivity of the catalyst and the selectivity of product s. Because these factors are closely associated with the surface processes such as adsorption, diffusion, and growth, as shown in Figure 1 2, a better understanding of these processes is imp ortant from an engineering perspective to govern the behaviors of a wide range of applications and from a fundamental science perspective to provide principle s for future design. 1.2 Integration of Computational and Experimental Methods The last two decades have witnessed a dramatic rise in computational resources that has facilitated tremendous progress in computational science. In particular, this progress has enabled the application of computational simulations to complement the experi mental observations for capturing th e process structure property re lationship of an interested material. V arious simulation methodologies have been developed to
16 elucidate materials behaviors over a wide length and time scale as shown in Figure 1 3. For ex ample, density fu nctional theory (DFT) which is a quantum mechanical modeling method, may be used to predict complex interfacial phenomena such as wetting, adsorption, diffusion, and segregati on for atoms, molecules and condensed phases. B ecause of its ex plicit treatment of electronic structure, i t is further able to provide high fidelity prediction s but unfortunately the computational cost increases rapidly with system size. Therefore, the se calculations are limited to a relatively small number of atoms (<500). To overcome this limitation, classical empirical methods (also known as interatomic potentials) that model materials at the atomic and nanometer scale without explicitly treating electrons have bee n developed and employed in molecular dynamics (MD) and Monte Carlo (MC) simulations. Such simulations have been employed to examine device sized systems and dynamic al process that include the effect s of temperatures and pressure, which is beyond the reach of quantum based approaches. The main strength of classical empirical potentials is their low computational cost relative to electronic structure calculations. Recently, computation of a system with over a trillion (10 12 ) atoms become s a reality in MD sim ulations 18 In MD simulations, the motion. Therefore, it is important that the interatomic potential, whic h is used to calculate energies and forces, is accurate and transferable in order to correctly describe the chemical and physics properties of materials under a variety of conditions. Over the last few decades, empirical potentials have been successfully d eveloped for the specific type of chemical bonding within a given system. For example Lennard Jones (LJ)
17 potentials 19 20 are typically used to characterize van der Waals (vdW) interactions ; embedded atom method (EAM) potentials 21 22 are typically used to describe metallic bonding; and Buckingham potentials 23 24 are typically used to model ionically bonded systems. Although these potentials and many others 25 26 describe the specific interactions for which they are parameterized ve ry well, they are typically applicable to materials with only one type of bonding which prevents their widespread application to a system with different types of bonding One example is heterogeneous catalysis as illustrated in Figure 1 4 which shows t hat catalytic driven chemical reactions involve severa l different interactions across interfaces such as reactant metal and metal metal oxide that may contribute to overall catalytic activity and selectivity. Therefore, a potential should no t only capture multiple bonding events and types but also correctly describe the ir relative strength. Recently, the charge optimized many body (COMB) 27 30 potentials and r eactive force field (ReaxFF) 27 30 potentials have been de veloped to tack le this challenge and parameterized to include a wide range of elements and compounds inclu ding heterogeneous systems The first generation COMB ( COMB1) potential was developed in 2007 by Yu et al 31 based on Tersoff potential 27 and the electronegativity equilibration (QEq) method proposed by Rappe and Goddard 32 33 Later, to improve the description of interfacial structures and s m all molecule s the second generation COMB (COMB2) potential in two different versions, COMB2A 34 and COMB2B 28 were developed. R ecently, the third generation COMB (C OMB3) 29 potential was developed usi ng terms from COMB2 potential 30 for electrostatic effects, the second generation reactive empi rical bond order
18 (REBO2) 35 potential for short range interactions, and the coordinat ion function that was developed for the REBO potential for MoS 2 36 systems In particular, because of the replacement of the original Tersoff expression, which lacks the four body dihedral term that is needed to capture the delocalized bonding in hydrocarbon systems, and with the short range interaction terms used in REBO2, the COMB3 potential can be extended to C/H/O/N systems 30 In this dissertation, the main objective is to apply and parameterize the COMB 3 potential for various heterogeneous systems and investigate the underlying physics of the interfaces and surface phenomena. The interfaces phenomen a will focus on the metal clusters growth on the metal oxide via deposition and the surface phenomenon will focus on the oxidation of metal nitride surfaces. 1. 3 Outline For a systematic study of this dissertation, the following chapters are or ganized as follows. Chapter 2 reviews the computational methods used throughout our investigated systems Specifically, the formalism s and parameterization procedures for the COMB3 potential are thoroughly described C hapter 3 presents the COMB3 potential development for the Cu/Zn and Cu/ZnO systems. Th is potential function is applied in MD simulations to model the deposition and subsequent evolution of Cu clusters on ZnO in order to elucidate aspects of the growth mode of Cu on ZnO that are currently the subject of controversy in the literature Chapter 4 provide s the COMB3 potential function for titanium (Ti) and titania ( TiO 2 ) The Ti/TiO 2 poten tial successfully describes the properties of hcp Ti metal and rutile TiO 2 In addition, the phase orders of TiO 2 polymorphs and surface energies of various rutile surfaces are
19 well predicted compared to DFT calculations The TiO 2 potential is further combined with the Cu/Cu 2 O and Cu/Ti COMB3 potentials to examine the energetics of different sized Cu clusters on rutile TiO 2 (110) surface. Chapter 5 provides the COMB3 potential function for TiN systems. The potential is further use d to study the structural and adhesive properties of Ti/TiN interfaces and characterize the adsorption behavior of oxygen atoms and molecules on the TiN (001) surface. Chapter 6 discusses behavior s of adsorption and subsequent dissociation of O 2 on the AlN stepped surface investigated by DFT calculation s. The results are further compared to behaviors of the flat surface. Lastly, the general conclusions of th is collection of works and future applications of the COMB 3 potential are given in Chapter 7
20 Table 1 1 Characteristics of physical adsorption and chemisorption. 13 Physical adsorption Chemisorption Weak van der Waals interactions S trong chemical bonding The heat of adsorption are low (20 40kJ/mol) The heat of adsorption are low (40 400kJ/mol) No electron transfer Electron transfer leading to bond formation No dissociation of adsorbed species May involve dissociation Nonspecific (insensitive to surface symmetry) High surface symmetry specific (ex: steps) Rapid, non activated, and reversible Slow, activated, and irreversible Only significant at relatively low temperatures Over a wide range of temperatures
21 Figure 1 1. Initial oxidation processes of bare metal surfaces.
22 Figure 1 2 Surface processes of metal thin film s or clusters growth on metal oxide substrate. 37 38
23 Figure 1 3. Hierarchy of multiscale modeling
24 Figure 1 4 Snapshot of a molecular dynamics si mulations of water, carbon mono xide, and carbon dioxide molecules interacting with a copper catalyst supported by a zinc oxide substrate. The element of Zn(grey), O(red), Cu(yellow), C(light blue) and H(white) is distinguished by color.
25 CHAPTER 2 COMPUTATIONAL METHODS 2.1 Computational Method s D ifferent types of computational methods exist that provide information regarding material properties from the electronic scale ( e.g., density of states) to the structural scale (e .g., stress concentration ). Ther e fore, understanding the advantage s and limita tions of different computational approaches is important for their proper application Here, the fundam entals of density functional theory (DFT) which considers interactions between atoms and electrons explicitly and classical atomic scale methods, which consider interactions between atoms with interatomic potential functions are discussed. These two app roaches are also the main tool s used in the computational work that forms the basis of this dissertation 2. 2 Density Functional Theory In density functional theory the many body electronic wave function of quantum mechanics is replaced with as the central qua n ti ty. The development of DFT is based on two theorems that were proposed by Hohenberg and Kohn in 1964 39 The external potential v ext ( r), and hence the total energy, is a unique functional of the electron density n expressed as a one to one relationship between the e xternal potential and the ground means that the energy of a m aterial has a functional dependence on single electron density, ( r ), and then ( r ) is a function of the positions of the electrons, r T he second The ground state energy cam be obtained variationally : the total energy is the exactly ground density words, the exact
26 ground state density of the system is that which minimize s the total energy. Therefore if we know the total energy functional in terms of the density, the ground state energy can be found by nu merically minimizing this function al. However, neither theorem defined the functional. Based on the Hohenberg and Kohn theorem s the electron density can be used to determine the ground state. However, we still need to determine the electron electron interaction term of the Schrdinger equation In 1965, Kohn and Sham 40 used the derivation of a set of equations that only involve single non interacting electrons. Therefore, the Schrdinger equation can be written in terms of the Kohn Sham equations as: (2 1) ; (2 2) The individual terms on the right side in Equation 2 1 are, in order, t he kinetic energy of electron, the electron nuclei potential the electron electron potential ( also called the Hartree potential), and the exchange correlation potential. is the wavefunction of state i that only depends on three spatial variable s and is the Kohn Sham eigenvalue. In E quation 2 1, the exact functional for the exchange correlation energy is the only unk nown term. Therefore, reasonable approximations such as the local density approximation (LDA) and generalized gradient approximation (GGA), are r equired to solve this equation. These approximations have been shown to work well for most
27 material systems, with the notable exception of strongly correlated systems 41 43 For more complete details regarding developments and application s of DFT, the reader is referred to Ref. 44 2. 3 Classical Empirical Method s 2. 3 .1 Interatomic Potentials In contrast to DFT, interatomic potential s are less computationally expens ive because they only consider interactions between atoms without explicitly treating electrons Through mathematical m odel s the energy of different types of bonding can be represented as function s of interatomic distance electronic density, bond angles, charge, types of neighbors, and/or other factors The most well known is the Lennard Jones (LJ) potential 19 20 which was initially proposed for liquid argon and takes the form of Equation 2 2: (2 3 ) where the r ij ( ) is the interatomic distance between two atoms and the parameters of and are associated with the materials property In addition the r 12 term represents the short range repulsion interaction from overlapping electron orbitals with equal spin (Pauli repulsion). In contrast, the r 6 term represents long range van der Waals attraction. A schematic curve of the LJ potentials energy ( U ) as a function of interatomic distance ( r ij ) is shown in Figure 2 1 (a). The parameter (eV) governs the strength of the interaction where the minimum potential energy min corresponds to an interatomic distance () In addition, the shape of the p otential energy curve controls many of the
28 physical properties of materials, such as the bulk modulus and thermal expansion. A schematic curve of th e force as function of interatomic distance is illustrated in Figure 2 1 (b), which are determined from the gradient of the interatomic potential energy and thus take s the form of Equation 2 3: (2 4 ) Form Figure 2 1 (b), it is clear that is the equilibrium interatomic distance where the net forces are zero. There fore, as indicated in Figure 2 1 (c), when the interatomic distance is larger or smaller than the equilibrium distance, an attraction or repulsion force acts and pushes both atom s until the y reach an equilibrium distance Because of its computational simplicity the LJ potential is extensively used In modern materials modeling, it is best applied to describe the properties of gase s or long range interactions Over the last few decades, empirical potentials using different mathematical expressions and various parameters for the potential function have been successfully developed for specific type s of chemica l bonding within a given system, as listed in Table 2 1 These potentials are then used to determine interatomic forces in for example, m olecular dynamics (MD) simulations. 2.3.2 Molecular D ynamics (MD) Simulations MD simulation is used to investigate the dynamical properties of a system at the atomic scale In this approach t he motion of each atom is described by solving Newt : (2 5 )
29 (2 6 ) where is the force exerted on the atom, which can be determined from the gradient of the interatomic potenti al energy In addition, is the acceleration of the atom and is its instan tan eous position Therefore both DFT and interatomic potential s may be used to determine interatomic forces within MD simulation s, although the latter is more computationally efficient and therefore more commonly used In addition, m any numerical algorithms have been developed for integrating the equation of motion, such as the Verlet, leap f rog, and predictor corrector. 45 In this dissertation, the Verlet algorithm is used. 2.3. 3 Ensembles MD is a statistical mechanics method, and it provides a way to obtain a set of configurations (atomic positions and momenta) distributed according to some statistical distribution function ( an ensemble). Therefore, in an MD simulation, various thermodynamic ensembles are employed. Three common e nsembles are: (1) Microcanonical Ensemble (NVE): used for an isolated system with N atoms, which keeps a constant volume (V) and a conserv ed total energy (E). (2) Canonical Ensemble (NVT): used for a system with N atoms in a temperature bath. The volume (V) and the temperature (T) of the system are kept constant. (3) Isobaric Isothermal Ensemble (NPT): used for an isolated system with N atom s in a temperature and pressure bath. The pressure (P) and the temperature (T) of the system are kept constant
30 2.3. 4 Periodic Boundary Conditions and Cutoff Radius At current levels of computing power, MD simulations are routinely carried out on system si zes on the order of several nanometers. However, they are useful in showing atomic scale details related to experimentally observed macroscopic systems. To mimic materials in the bulk or at surfaces, periodic boundary conditions (PBCs) are typically applie d. In particular, three dimensional PBCs are used to mimic a p ortion of a bulk system with no surfaces present. In contrast, simulations of planar surfaces make use of two dimension al PBCs within the plane of the surface ; the third direction typically incl udes a vacuum region. In general, large systems can be more readily modeled with smaller atomic scale systems through the use of PBCs. Figure 2 2 illustrates PBCs in two dimensions. The system cell is repeated to mimic an infinite lattice and thus does not have edge effects During the simulation when a particle moves in the system cell, the corresponding particle in each neighboring periodic cell moves in the same manner. If a particle leaves the system cell, say through the left boundary, one of its perio dic images enters throu gh the opposite side, the right boundary. 2. 4 Variable Charge Reactive Potentials: COMB Potential s 2. 4.1 General Formalism of the Third G eneration COMB P otentials This section presents the general formalism of COMB potentials and describes how they are used to obtain a proper description of energy and charge for a multicomponent system As show in E q uation 2 7 the total potential energy ( ) of a system as described by the COMB potential is composed of electrostatic energy ( ) short range interactions ( ), van der
31 Waals interactions ( ), and correction terms ( ), where q and r represent the charge and coordinate array of the system, respectively. (2 7 ) 2.4.2 Electrostatic Energies The electrostatic energies ( ) include t he self energy ( ), the charge charge interactions ( ), the charge nuclear interactions ( ), and the energies associated with atomic polarizability ( ), which are defined in Equation 2 8 The repulsion between nuclei is not explicitly formulated since their overall effects ar e assumed to be included in the short range interactions (2 8 ) The self energy ( ) descr ibed the energy required to form a charge on an atom. This can be expressed by the summation of the ionization or electron affinity energy of an isolated atom ( ) and a correction function ( ), termed the field effect 46 that reflects the change of electronegativity and atomic hardness of the atom with its e nvironment, as specific in Equations 2 9, 2 10, and 2 11 : (2 9 ) (2 10 ) (2 11 )
32 In Equation 2 10 as outlined by Toufar et al. 47 the parameters is identified with the electronegativity and J is associated with chemical hardness or self Coulombic interaction 47 which are inherent to each element. Therefore, i J i K i L i are treated as atomic parameters and fitted to the ionization energies and electron affinities for each of the elements. In contrast, the parameters and in Equation 2 11 are associated with specific bonds to describe the environmental dependence of field effects. The long range electr ostatic interaction between two ions is normally described which is defined in Equation 2 12 (2 12 ) However, as the distance between two ions approaches to zero, the Coulomb expression in Equation (2 6) results in infinite values, which this Coulomb catastrophe that was introduced by a point charge model, COMB potentials adopt the Streitz Mintmire 34 charge density function. Specifically, the charge density of an atom is taken to be function of its charge, spatial location ( r ) and atomic position ( r i ): (2 13) (2 14 ) is an effective point core charge treated as a fitted parameters, is the Dirac delta function, is a function that capture the radial decay of the electron density of the S type orbital 48 and is an orbital exponent that controls the length scale associated with this decay. Further, the charge charge interactions ( ) can be calcu lated as
33 indicated in Equation 2 15 by the integration over electron densities between pairs of atoms through the Coulomb integral, defined in Equation 2 16 (2 15 ) (2 16 ) W here and indicate the centers of and and is the distance between density distribution. Similarly, the charge nuclear interaction are expressed in Equation (2 1 7 ) through the charge nuclear coupling operator, which is shielded nuclear attraction in tegral and defined in Equations 2 18 and 2 19 (2 1 7 ) (2 1 8 ) (2 1 9 ) The Coulombic intera ction in Equation (2 9) and (2 11) are solved via a direct Wolf summation method 49 which is more computationally efficient than the Ewald summation 50 for computing electrostatic interactions. The polarizability, P presents the relative tendency of an equilibrium charge distribution to be distorted in response to an extern al field, which is defined as the ration of the induced dipole moment, of an atom to the electric field, The inclusion of the polarization effects has been proven to be useful in classical simulations to stabilize alumina 51 and to model the interactions with
34 small polarizable molecules such as water and oxygen molecules. Using the fluctuation dipole model by Sterne et al. 52 as specific in Equation 2 20 the polarization vector can be directly calculated from the ele ctrostatic field generated by the atomic charges, and the neighboring induced diploes. (2 20 ) (2 21 ) is the polarizability tensor and is the induced dipole field tensor. Because the atomic polarizability is assumed to be isotropic, thus reduces to a scalar value that is only associated with the elemental nature of the atom. The induced dipole field tensor ( ) is employed as a damped function that diminishes as atoms overlap. (2 22 ) Lastly, as shown in Equation 2 23 the energies contributed from atomic polarizability are the dipole self energy, the dipole charge interaction, and the dipole dipole interaction, respectively. (2 23 ) 2.4.3 Short R ange interactions The bond order type short range interactions, where the short range repulsion energy, and attraction energy, are based on the
35 Tersoff 53 and Yasukawa 54 po tentials, as briefly described in Equation 2 24 For the origin Tersoff potential, the short range interaction is dependent on the inter atomic distance ( ). However, the short range interaction in COMB potentials is not only exponentially decay with interatomic distance ( ) but vary with charge ( ). In Equations 2 25 and 2 26 the charge dependent correction functions, are added to the exponential t erm of the repulsion and attraction energies to reflect the change in atomic radius with charge. The change in the bond order with charge is re flected in the charge dependent function, In addition, the cut off func tion, as descr ibed in Equation 2 27 adopts the cosine decay formalism to smoothly terminate the short range interactions, where and are the upper and lower cutoff distance. (2 24 ) (2 25 ) (2 26 ) (2 27 ) The bond order term ( ) uses a similar formalism that is used in the REBO potential 55 and includes the contributions from the bond angle ( ), coordination ( ) additional torsion ( ), and conjugation ( ) effects from hydrocarbon
36 systems, which are used to capture many body effects. Here, due to the complex bond environment of hydrocarbon system, th e bond order term only show the contri bution from bond angle and coordination without torsion and conjugation effects, as shown in Equation 2 28 The detailed parameters and fitting procedure of the COMB potential for hydrocarbon system is described in ref 36 56 (2 2 8 ) The asymmetric term is used to weaken the longer bond. The angular function, and the coordination function, are described from Equations 2 29 to 2 31. (2 2 9 ) (2 30 ) (2 31 ) In Equation s 2 2 9 to 2 31 to and to are fitted parameters, and is the number of nearest neighbors of atom 2.4.4 v an d er Waals Interaction s The long range van der Waals (vdW) interaction s ( ), are captured by the classical Lennard Jones formula in COMB potential, and take the form of Equation 2 32 : (2 32 )
37 In this expression, and reflect the strength and equilibrium distance of the vdW interactions, respectively. The vdW interactions are truncated and shifted to zero at the Coulombic cut off radii. To avoid extremely high repulsion between short range bonded atoms, a cubic splin e function is added to smoothly terminate the vdW interactions at the upper end of short range cut off radii. For any binary vdW bonds, the and are decided by considering the geometric and arithmetic m eans of values of element type vdW bonds, respectively. 2.4.5 Correction T erms The correction terms, ( ) are primarily used to modify the energy contribution from specific bond angles. In the set of correct t erms, as described in Equation 2 33 they consist of Legendre polynomials ( LPs ) up to sixth order and a bond bending ( BB ) term. (2 33 ) The LPs whose detailed forms can be found in ref. 30 are a set of symmetric energy penalties on the bond angle. In contrast, the bond bending term pro vides an asymmetric energy penalty on a specific bond angle, For more detailed on the application of correction terms in COMB potentials, they can be found in refs 57 All the parameters in Equation (2 27) are thre e dimensional, where the first subscript character ( ) represents the element type of the central atom. 2.4.6 Dynamic Charge Equilibration Method The key development of COMB potential applicable to multiple types of boding is a met hod for equilibration of charge. Although fixed charged potentials work well for
38 bulk materials systems, they cannot describe systems that require charge equilibration, such as a heterogeneous interface between dissimilar materials. In COMB potential, the atomic charges are treated as dynamical variables evolving with time and are equilibrated based on the electronegativity equilibration (EE) principle. This principle, which was proposed by Sanderson, states that in a closed system of interacting ions at ch emical equilibrium, the electron density is distributed so that electrochemical potential is equal at all atomic sites. Base on the definition that the chemical potential ( ) for each site (atom) is equivalent to the partial deriva tiv e of the total energy ( ) with respect to electron density ( ) : it is also equivalent to the negative of the electronegativity ( ) at each site. Therefore t he partial charge on each atom can be determined as that which corresponds to the equilibrium electronegativity To efficiently and explicitly determine the pa rtial charge and position of each atom dependent upon time, an extended Lagrangian me thod developed by Rick et al. is employed in COMB potential. The Lagrangian for the system is: (2 34 ) In Equation 2 34 and correspond to the mass, position, and charge of atom respectively. The first term on the right side is the kinetic energy associated with the physical movement of the ions. The second term is the kinetic energy of the charges with a fictitious charge mass ( ). The potential energy is defined in third term. The final term is the contrast, where is an undetermined mu ltiplier that enforces the const raint to keep the charge in the system conserved.
39 From the Lagrangian, we can derive equations of motion for the atomic positions and charges as: (2 35 ) (2 36 ) Since the total charged is conserved in a system, it follows that: (2 37 ) Therefore, solving the above two equations shows is equal to the mean electronegativity for th e system, as shown in Equation 2 38 ( 2 38 ) This leads to Equation 2 36 can be rewired as Equation 2 39 : (2 39 ) The partial charge on each atom is determined by iteratively solving the equations of motion with damping factor for each set of nuclear pos itions until the electronegativity ( ) on each ( ) of the system. 2.5 Parameterization of COMB Potentials Although the COMB potential has a comp lex formalism for the potential energy, each component is indispensable to model multiple types of bonding within a given system. Based on the origin of the potential energy, the general procedure for parameterizing the COMB potential is performed as follo ws 30 58 : (1) creating the
40 database from published experimental data and/or electronic structure calculations; (2) parameterizing the pure systems ( starting with the parameters of short range interactions and the n that of electrostatic energies); (3) parameterizing the binary systems (starting with the parameters of short range interactions). 2.5.1 Creating the Database The parameterization of COMB potential proceeds by optimizing its parameters against the fittin g database. In order to allow it to be suitably transfer rable between various systems, a fitting database that includ es a wide ranges of pure and binary compounds of various phases needs to be considered in the parameterization process. Normally, the fitti ng database for pure system consists of bond length (or lattice constant), elastic c onstants, cohesive energy and for binary systems consist of atomic charge and heat of formation. The fitting database should also comp rise data from systems with a variety of coordination numbers to provide a good description of the bond order dependence, which therefore consists of different defected and surface structures. In addition, for different focus of studies, the database such as transition state energy for atomic diffusion can be added into parameterization process. The properties of these systems are obtained from published experimental data and/or the electronic structure calculations. For crystalline structure, the electronic structure calculations are performed with plane wave DFT calculation using the software of Vienna ab initio Simulation Program (VASP) 59 62 with appropriate pseudopotential and exchange correlation function (US LDA or PAW PBE). For molecular systems, the electronic struct ure calculations are performed using the Gaussian09 63 computational chemistry software package.
41 2.5.2 Parameterization of Pure Systems Because the electrostatic energy terms for an uncharged pure system (such as bulk strictures) is zero, the energy contribution of COMB potential for such a system is only from the short range interactions (Equation 2 18) and the form alism is reduced to the Tersoff ty pe of potential. In this case, this is straightforward for the parameterization of COMB potential, which only fitting parameters in the pairwise term and bond order term, as shown in Table 2 1. The van der Waals interaction will be considered in hydrocarbo n systems. After fitting the parameters in the short range interactions, the second step is to determine the charge associated parameters of an atom. The involved parameters for this step are listed in Table 2 2. The fitting process for this step is achieved through the following three sub steps: (1) fitting to the electron affinity and first, second, and third order of ionization energies of an isolate atom; (2) fitting the atomic polarizability, from the atomic polarization in a bonded dimer system; (3) fitting the rest of the parameters to the reaction energies of inserting or detaching an electron form the dimer sy stem 2.5.3 Parameterization of Binary Systems Lastly, the COMB potential is fit to the binary system with multiple phases. Similar to the parameterization process for pure system, the potential starts with fitting the pairwise terms to the phase orders and charged of variety o f binary systems, which is followed by fitting the parameters in many body terms for the properties of the phase of interest. The involved parameters for the binary systems are listed in Table 2 3.
42 2. 6 Transition State Search: Dimer Method Long time simulations of system evolution require approaches other than traditional MD, such as a kinetic Monte Carlo ( K MC) scheme. In a typical K MC simulation, the energetics of all transitions that might occur in the system must be known. By contr ast, methods such as the s tring approach 64 and nudged elastic band (NEB) 65 method are used to identify transition states. However, the requirement of knowing the initial and final state limits the application of these approaches to simple systems. In the present work, to further evaluate the capability of the COMB potential and to characterize the transition states of Cu motion on the ZnO surface, the d imer method developing by Henkelman and Jnsson 66 is utilized. In the dimer method algorithm, effective forces are used to move along the lowest curvature dire ction in order to find the saddle point. Because the dimer method uses only the first derivatives of the potential and the initial state, it is computationally efficient. The diffusion Cu adatom on the Cu(100) surface has been investigated extensively in experiments 67 68 and simulations 69 72 Therefore, this system is used as a first test of our implementa tion of the dimer method with the COMB potential. In MD simulations with the embedded atom method (EAM) potential, Liu et al. 70 predicted the migration barriers for the hopping and exchange mechanisms of Cu on Cu(100) to be 0.38 eV and 0.72 eV, respectively. Hanson et al. 69 also examined these mechanisms with DFT calc ulations and found that the hoping mechanism had an activation energy of 0.69 eV, compared to 0.97 eV for the exchange mechanism. In the present work the Cu (100) surface consists of 864 atoms, 72 atoms per layer in 12 layers. A single Cu adatom is placed in 4 fold hollow site which is the site with the highest adsorption energy. Approximately 20 of vacuum is used to prevent interaction between the top
43 layer and a periodic image of the bottom layer. In addition, only the Cu adatom and first layer atoms a re allowed to move freely. Figure 2 3 shows the two transition mechanisms for hopping and exchange as predicated by the COMB potential. The predicted migration barriers are 0.42 eV and 0.53 eV, respectively, which are in qualitative agreement with experime ntal findings, prior EAM simulations and DFT calculations. This test illustrate that the COMB potential is able to predict the correct transition mechanisms of Cu on the Cu (100) surface. 2. 7 Summary Generally, the COMB potential has thirty elem ent typed (or called one dimensional) parameters for each element, thirty two bond typed (two dimensional) parameters for each bond, and additional bond angle typed (three dimensional) parameters, which is optional, for correcting the energy contribution from speci fic bond angles. This design is important to allow parameters to be flexible and transferable in a multicomponent system. Therefore, for the COMB potential, there is no re parameterization needed of previous parameters when the potential is expanded to inc lude more components or when new binary systems are examined Therefore, relatively speaking the parameterization of COMB potential s is more challenging than traditional interatomic potentials because of its considerable fitting database and parameters. W e have therefore developed an in house code called Potential Optimization Software for Materials (POSMAT) 58 to speed up the parameterization process by utilizing a shared memory parallelization scheme. In addition the COMB potential is currently implemented in the Large scale Atomic/Molecular Massive Parallel Simulator (LAMMPS) 73 packages distributed by Sandia National Laboratories. This availability allows the COMB potentials being capable of simulating system larger than
44 1 0 6 atoms at nanosecond timescale. A comparison of the performance of the COMB potential and other many body potentials within LAMMPS can be found in R ef 74
45 Table 2 1. The common used interatomic potentials. 75 Potentials Formalisms Application Lennard Jones 20 (1924) The constants of and are dependent on the physical property of materials. Mostly suitable for rare gases Buckingham 23 (1938) The constants of and are dependent on the physical property of materials. Mostly suitable for covalent bonding Embedded Atom Method ( EAM ) 22 (1984) is the total electron density at the atom due to the res t of the atom in the system. are constants dependent on the materials. Mostly suitable for metals Tersoff 32 (1988) and are the potentials due to repulsive and attractive forces between atoms and ; is a bond order parameter which takes the many body effects into account. Mostly suitable for covalently bonded materials R eactive empirical bond order (REBO) potential 56 ( 1990) Mostly suitable for covalently bonded materials and hydrocarbon
46 Table 2 2 List of fitting parameters for the pure systems. Parameters Descriptions Equation s pairwise term 2 19 2 20 in bond order term 2 28, 2 29 asymmetric term 2 29 to for bond bond angle term 2 30 to for bond coordination term 2 31 vdW interaction 2 32 to and correction functions 2 33 Table 2 3 List of fitting parameters for the charge dependent energy terms. Parameters Descriptions Equation ionization energy 2 1 0 field effects 2 11 Coulomb interaction 2 13,2 14 atomic polarizability 2 20 charge dependence on short range 2 25, 2 26 Table 2 4 List of fitting parameters for the binary system s. Parameters Descriptions Equation field effects 2 1 1 pairwise short range 2 19,2 20 symmetric term 2 20 to for and bonds bond angle term 2 30 to for and bonds coordination term 2 31 to and correction functions 2 33
47 Figure 2 1 Lennard Jones (LJ) potential. A) Energy function B) f orce function C) illustration of repulsion and attraction forces when an atomic distance is below or over the equilibrium distance.
48 Figure 2 2 Schematic representation of pe riodic boundary conditions in MD simulations. The square with dark background represents the simulation cell, and the others are the images.
49 A B Figure 2 3 The two migration mechanism s of single Cu atom on Cu(1 00 ) A ) H op B) e xchange. The top schematic views, shown from left to right, are the initial state, saddle points, and final state, respectively. The yellow spheres are sublayer atoms, the orange spheres are surface atoms, and the green sphere is an adatom.
50 CHAPTER 3 DEVELOPMENT OF COMB3 CU/ZNO POTENTIAL AND ITS APPLICATION ON CU CLUSTERS DEPOSITION ON ZNO SURFACE Methanol (CH 3 OH) is a potential candidate for the next generation of renewable green fuels because its energy density is sufficiently close to that of gasoline. In addition, methanol exists in liquid form at ambient pressures and temperatures, which makes it easier to store and transport than hydrogen. Methanol synthesis from the hydrogenation of CO 2 is performed commercially using Cu/ZnO/Al 2 O 3 catalyst. Despite numerous experimental and theoretical studies, there are still questions about the nature of the active site of the Cu/ZnO catalysts. Cu metal clusters 76 Cu + 77 80 the Cu/ZnO interface 81 83 and Cu Zn 84 alloys have all been proposed as the active site. In addition, different sites might be active for CO 2 hydrogenation versus CO 82 83 85 While Cu/ZnO has been used in syngas hydrogenation, there is recent interest in using such catalysts in electrochemical reduction of CO 2 to produce methanol. This interest is motiva ted by studies that show that Cu oxide can electrochemically reduce CO 2 to methanol 85 Le and co workers have correlated methanol formation on the Cu oxide elec trodes to the presence of Cu + sites 8 86 87 However, Cu oxide electrodes are inherently unstable in the reducing condition and methanol producti on tapers off with time as the oxide is reduced to metal. Recent results suggest that Cu on ZnO could lead to higher methanol formation rates but the nature of the active site and the stability of such catalysts still need to be explored 8 Since the morphologies of Cu cluster on ZnO surfaces can undergo dynamic changes, it is chall enging to experimentally specify the state of the active sites under reaction conditions. Thus, to elucidate the formation of methanol on Cu/ZnO catalysts, the first step is to characterize the behavior of Cu atoms on ZnO surfaces. There have
51 been several surface science studies to characterize Cu growth and stability on polar 88 and nonpolar ZnO surfaces 37 89 93 The polar ZnO surfaces are known to undergo complex reconstructions 94 and are therefore more difficult to examine using density functional the ory (DFT). Consequently, the majority of the DFT work 95 has been focused on understanding Cu behavior on the ZnO and ZnO surfaces. Hu et al. 96 97 investigated Cu deposition and subsequent growth on the ZnO and ZnO surfaces using DFT. Their results indicate that Cu atoms strongly interact with ZnO becoming slightly positively with a charge that between those suggested by Klier 97 (Cu + ) and by Fleisch 98 (Cu 0 ). It was suggested that the charged Cu atoms and the Cu Zn sites might be the active sites of the system. While first principles methods, such as DFT, provide high fidelity because of their explicit treatment of electronic structure, they are limited in the size and time scales that can be addressed directly. Because of the potential importance of lar ge scale changes of the Cu/ZnO interface (sintering of Cu metal, alloying with the oxide surface, oxidation of Cu metal, and adsorbate induced surface reconstruction) on the catalytic activity, a computationally more efficient method is needed. One potenti al option is the development of potentials (reactive force fields) that are able to capture charge transfer within a system (such as at the Cu and ZnO at the interface) but still allow for much more rapid computation than DFT methods. Recently, the reactiv e force field (ReaxFF) method 86 has been applied to model related complex reactive systems including ZnO/water 99 and Cu/Cu 2 O/water 100 101 However, as yet these potentials have not been directly applied to the Cu/ZnO system.
52 Here, a third generation charge optimized many body (COMB) potential for ZnO is developed and used to model the ZnO surface interacting with either a single Cu adatom or a beam of deposited and equilibrated Cu clusters in classical molecular dynamics (MD) simulations. The use of a reactive potential with dynamic charge 102 such as COMB, in the simulations under temperature and pressure conditions that are similar to those in the experiments is necessary to fully describe the potential nanometer scale changes at the interface, including sintering of the Cu clusters, al loying with the oxide surface, and oxidation of the Cu during the growth process. The COMB predictions for ZnO properties and of adsorption and migration energies of Cu atoms on the ZnO surface are compared to published experiment al data and to the results of density functional theory (DFT) calculations. In addition, the MD simulation results are compared with STM results to explain some puzzling experimental observations. The simulations further explore the way in which the growth mode and final morphology of the Cu cluster s are influenced by cha n ges in the Cu incident deposition energy and surface temperature. Cu growth on ZnO is also compared to Cu growth on Cu to determine how the nature of the surface influences the predictions 3 .1 Potential development and Properties of Cu/ZnO 3 .1.1 Fitting Procedure To ensure the phase stability of ZnO polymorphs, in addition to the wurtzite phase, another three phases of ZnO were explored: zincblende ( CN Zn O : 4), rocksalt ( CN Zn O : 6), and caesium chloride ( CN Zn O : 8) were considered. Since the bond lengths (Zn O: 1.99 vs 2.00 ) and bond angles (Zn O Zn: 108.55 0 and 110.38 0 vs 109.47 0 ) of wurtzite and zincblende are similar, there is only a small difference in cohesive
53 energy between tw o phases. The DFT PBE calculations predict that the energy of the zincblende phase is only 0.014 eV per ZnO unit higher than that of wurtzite phase, which is sim ilar to the energy difference calculated using the local density approximation (LDA) and genera lized gradient approximation (GGA) (about 0.015 and 0.013 eV, respectively). The parameters for ZnO COMB potential were obtained through fitting to the data set including the energy dependence on volume of four ZnO polymorphs, elastic constants, and surfac e energies from experiments and DFT calculations 3 .1 .2 Predicted Properties of ZnO polymorphs Table 3 1 compares the properties of lattice constants, elastic constants, and cohesive energy for various ZnO polymorphs predicted by the COMB2 and COMB3 potentials with values from experiments, DFT calculations and other empirical potentials. I t is clearly shown that the newly developed COMB3 potential exhibits good predictions as COMB2 potential on the cohesive energy and the elastic constants of the wurt zite phase. Most importantly, the incorrect surface energy predictions of the COMB2 potential have been improved in the new COMB3 potential, where the polar surfaces of wurtzite ZnO have higher cleavage energies than nonpolar surfaces. Nevertheless, becaus e the stabilization effect for the ideal termination of the polar surfaces of ZnO involves charge compensation of the surface layers, the cleavage energy for the polar surfaces is not quantitatively captured in the current COMB potential. Therefore, in fut ure work incorporating data related to the polar surfaces, such as cleavage energies or reconstructed surfaces, may improve the COMB potential description of the polar surfaces.
54 3 1.3 Predicted Properties of CuZn alloys The Cu Zn interactions in COMB 3 are explicitly parameterized for three specific Cu Zn alloys 103 104 which are the brass Cu 3 Zn (fcc), brass CuZn (bcc), and brass Cu 5 Z n 8 ( cI 52) phases, respectively. Table 3 2 lists the COMB 3 predations and makes a comparison to DFT and published experimental data. The heat of formation for the phase is properly calculated to be the lowest and the phase order prediction is correct. experimentally, the correct phase order is of more importance because of the foc us on potential alloying during cluster deposition and metal thin film growth rather than mech anical properties of this alloy. 3 2 Interaction of Cu on ZnO 3 .2 .1 Adsorption Energies of Cu Adatom on ZnO Several computational studies 105 of Cu interacting with the ZnO surface have been carried out. For example, Beltrn et al 96 97 106 107 examined Cu adsorption directly on top of O and Zn surface atoms using Hartree Fock theory and predicted that adsorption on O was slightly more favored than adsorption on Zn. Subsequently, Hu et al 96 and Hellstrm et al 97 performed more extensive DFT calculations that allowed for geometry optimization and predicted three poss ible adsorption sites for Cu on ZnO as illustrated in Figure 3 1 The ir calculation s further predict ed that site 2 in Figure 3 1 was considerably more stable than sites 1 and 3. These three adsorption s ites are investigated here with the COMB3 potential. Table 3 lists the adsorption energ ies and charge s of single Cu atom s at the s e sites ; a more positive value corresponds to a more stable configuration. The COMB 3
55 potential predicts that the most favorable adsorption site is site 2 and the least favorable is site 1 which is consistent with the published DFT results 107 at low Cu coverage. In addition the atomic charge s on the absorbed Cu is predicted to be higher at site 2 than at the other two sites which correlates with the strong interaction between Cu and two O atoms. 3 .2 .2 Migration barriers of Cu Adatom on ZnO To characterize the migration of a single Cu atom on the ZnO surface, the dimer method 107 was employed to characterize its transition path. Table 4 lists the barriers for Cu atom migrat ion from site 2 to sites 1 and 3 Although the barriers predicted by COMB 3 are higher than those predicted by DFT 108 the two sets of barriers are in qualitative agreement with one another. Specifically, both predict that Cu atom migration along the direction is more energetically favorable than along the direction. The relative stability of different sized copper clusters on the ZnO surface can be characterized by the binding energy. The binding energy ( E b ) per atom may be calculated using the following standard expression, ( 3 1) where n is the number of Cu atoms, E ( Cu n /ZnO ) is the total energy of Cu n adsorbed on the surface, and E ( Cu n ) and E ( ZnO ) are energies of Cu n clusters and the optimized bare ZnO surface, respectively. The more positive value of E b correspond s to a more stable species. Figure 3 2 illustrates the binding energy for Cu n (n=2 5) clusters predicted by DFT calculations and COMB potenti al, which shows a fairly good correlation between
56 the values from DFT and COMB based on the estimation of linear trendline (R squared value is 0.983). This result also indicate s that not surprisingly, Cu prefers to exist as a cluster than remain as an in dividual atom on this Zn O surface. Figure 3 3 shows the charge configuration of flat and planar Cu 5 clusters on ZnO which indicates that the charge dependence of Cu clusters in different configurations. From the above analyses it is clear that the COMB potential correctly describes the adsorption and migration of Cu atoms on ZnO Therefore, in the next section, this potential is used to elucidate Cu thin film growth via depositio n for homogeneous and heterogeneous interfacial systems. 3 .3 Cu Cluster s Deposition on ZnO 3 .3.1 Details of the MD Simulation of Cluster Deposition and Relaxation Classical MD simulations are used to model Cu 6 deposition on the ZnO surface, where the forces on the atoms are determined via the COMB 3 potential. Figure 3 4 shows the deposition system which is composed of deposited Cu cluster s and a ZnO surface slab. The dimensions of the ZnO surface slab are 72 74 70 and it consists of 36960 atoms; the direction is perpendicular to the plane of the surface. Periodic boundary conditions are imposed in the and directions to simulate an infinite surface. There are three types of atoms in the slab: fixed, thermostat, and active atoms. The fixed atoms comprise the bottom four layers (~ 10 ) and are held fixed in their bulk lattice sites to prevent the surface slab from translating as a result of the cluster deposition process. The middle fourteen layers (~ 30 ) of ZnO are thermostat atoms, and have a Langevin 107 thermostat applied to them to dissipate the excess energy from the deposited clusters and maintain a constant surface
57 temperature of 300 K. The active region comprises the top eight layers of th e slab and evolves under unc onstrained equations of motion. equations of motion, which are solved by using a Verlet algorithm 109 with a time step of 0.2 0 fs. The incident beam consists of Cu 6 with a stable octahedral structure that are deposited at a rate of 1 5.1 atoms/ps which is on the same order as the experimental value of 10.3 atoms/ps 110 The kinetic energy of the clusters is 1.25 eV per atom, and the fluence is 4.5 J /cm 2 All of the MD simulations are performed using the LAMMPS 111 software 3 .3 2 Surface Morphology and Growth Mode Figure 3 5 (a) i llustrates the morphology of the ZnO surface following the deposition of the equivalent of one monolayer (ML) of Cu which corresponds to 600 Cu atoms. The surface Cu atoms are colored by their height above the ZnO surface to distinguish the monolayer high 2D and multilayer high 3D Cu cluster s from one another as illustrated in Figure 3 5 (b). The MD simulations clearly predict that both 2D and 3D Cu cluster s are formed following deposition. To clarify the mechanisms associated with thin film growth, Figur e 3 6 illustrates the way in which the surface morphology of the Cu varies with coverage from 0.2 to 1 ML. Initially, at 0.2 ML, the deposited Cu atoms mostly maintain their cluster form on the surface while sometimes break ing into small er clusters or even individual atoms following their initial impact on deposition. Because the incident energy is able to provide the energy that adatoms require to overcome migration barrier s on ZnO
58 adjacent Cu adatoms are able to migrate along the direction and further coalesce into larger 2D Cu cluster s. When the Cu coverage reaches ~0.4 ML, a substantial fraction of the surface is cov ered by 2D Cu cluster s. However, the formation of a second layer of Cu is also pred icted, as ind icated by the circle in Figure 3 6 ; this is the start of 3D structure formation At this stage, the simulations predict that the deposited Cu clust ers increase the area of the 2D clusters when they land on either the edge of existing 2D cluste r s and diffuse down to the lower terrace, or on the clean ZnO surface as indicated in Figure 3 7(A ). In contrast, when the deposited Cu land around the center of existing 2D clusters they stick. This leads to the formation of the second laye r Cu atoms, a s shown in Figure 3 7(B ). These two behaviors may be further explained by the competition between diffusion energy and the Ehrlich Schwoebel ( ES ) barrier 74 76 The ES barrier, which is illustrated in Figure 3 8 E ES ) t hat is added to the diffusion barrier ( E diff ) as a result of bond stretching at the step when an adatom cross es the step edge to the lower terrace It is therefore easier for Cu atoms at the step edge to reach the lower terrace while, because of the longer path to the edge, Cu atoms in the ES and E diff and so are reflected from the step edge to form a second layer. As the Cu coverage increases above 0.6 ML, the spreading 2D clusters are found to exhibit a Cu (111) like orientation, as indicated by the triangle in Figure 3 6 As the Cu coverage increases, the 2D clusters grow increasingly thicker and more three dimensional with out filling in the gaps between the 2D clusters, such that the
59 percentage of uncovered ZnO surface remains essentially unchanged This implies that at this coverage, the 3D growth mode is dominant for Cu The growth of different metals on oxide surface s has been studied and summarized in several review articles 112 113 Growth mode s similar to those described here w ere first observed by Ernst et al 4 5 114 for Cu particles that were vapor deposited on oxygen terminated ZnO(0001) at 130 K. Th is study also found that the vapor deposited metal gr e w in the form of 2D clusters up to a certain critical coverage, beyond which 3D clusters form ed This type of growth mode was su bsequently observed in other Cu/ZnO systems 115 However, because the temperatures of the ZnO surface are below 300 K in the current work the growth mod e and surface morphology are likely kinetically limited Therefore, in the following sections, the kinetic effects are examined separately with respect to incident energy and the temperature of the ZnO surface. 3 .3 3 Effect of Incident Energy The incident energy of the Cu clusters is one of the parameters varied for the deposition process. A previous experimental study 4 92 115 117 of Cu thin films grown on Cu(111) by direct ion deposition exhibited a strong morphological dependence on the incident energy. In an effort to p robe this effect for Cu growth on ZnO, three different incident energies, 0.30, 1.25, and 2.81 eV/atom, are considere d for C u 6 up to a coverage of 1 ML while the surface temperature is maintained at 300 K. The resulting Cu thin film s are illustrated in Fig ure 3 9 At energies of 0.30 eV/atom, the Cu 6 octahedron clusters are predicted to spread out on the surface to form planar clusters or break into individual atoms without maintaining their origin al geometry. Because the incident energy is sufficient to overcome the migration barrier, adjacent Cu atoms or small clusters are able to merge
60 into larger clusters. This spreading and coalescence is predicted to be more prevalent at energ ies of 1.25 eV/atom which lead s to a larger percentage of the ZnO su rface covered by 2D Cu clusters. As the incident energy is increas ed further to 2.81 eV/atom, the incident Cu clusters are predicted to penetrate into the su rface and to partially mix with the surface Zn at oms as they rebound, thus forming Cu Zn mixtures at or near the surface. In an attempt to quantitatively analyze the se modifications to the surface morphology, the surface roughness of the deposited film is investigated by the approach of root mean square (RMS) method The definition of RMS can be expressed as (3 2) where Z i represents the height of the deposited atoms on the su rface denotes the mean height of all the deposited atoms, and n is the total number of deposited atoms. The surface roughness based on RMS is 1.48, 1.33, and 1.42 for 0.30, 1.25 and 2.81 eV/atom, respectively. This finding indica tes that the surface morphology becomes smoother as the incident energy increases up to a certain point. However, at energies beyond this point, the surface atoms will be sputtered and the surface roughness will increase. 3 .3 4 Effect of Surface Temperature In addition to the incident energy, the surface temperature is another deposition condition that can be varied. We considered the deposition of 1ML of Cu with an incident energy of 1.25 eV/atom on surfaces maintained at tempe ratures of 300, 600 and 800 K. As shown in Figure 3 10 although these predicted surface morphologies
61 finally exhibit 3D Cu cluster formation, small and random 3D clusters are found when the surface temperature is maintained at 800 K. Based on our simulati on results, the area of 2D Cu clusters in terms of the number of Cu atom in the first layer increases with increasing surface temperature. It is estimated that the area of 2D Cu clusters increases around 12 % when the surface temperature is 800 K relative to when it is 300 K. Additionally, the structures of the 2D Cu clusters are more compact and exhibit the well defined, Cu (111) structure at 800 K relative to their structures at 300 K. This relative dependence of surface temperature and Cu cluster morphol ogy can be attribute d to the fact that the thermal annealing effect provides additional energy for the migration of the Cu adatoms. However, because the thermal energy (0.08 eV at 800 K) is relatively small compared to the migration barrier (0.7 eV), the e ffect of surface temperature on the resulting surface morphology is not as significant as the effect of incident energy. 3 .4 Comparison of the Growth of Cu on Cu(111) and ZnO It is expected that the surface structure of the substrate will influence the behavior of the deposited Cu clusters. To fully explore this effect on Cu thin film growth, the parameters used for 2ML Cu deposition on the ZnO surface are applied to Cu deposition on the Cu(111) flat surface. Figure 3 11 and 3 12 illustrates the evolutions of surface morphologies on both surfaces. Based on our previous study 38 COMB predicts that the two lowest migration barriers for hopping and exchange mechanisms of a Cu adatom on Cu(100) are 0.42 and 0.53 eV, respectively. Because Cu (111) is a close packed surface the migration barrier for a Cu adatom on the (111) surface is lower than on the (100) surface. Thus,
62 as expected, our simulation results predict that deposited Cu atoms are able to easily migrate on the Cu(111) surface through hopping and exchange mechanisms which leads to the formation of laminar thin films with low surface roughness, as shown in Figure 3 13 (A). The result is in excellent agreement with experimental findings 106 In contrast, as show n in Figure 3 13(B) deposited Cu clusters are predicted to form larger 3D clusters companied by spreading 2D clusters on ZnO that exhibit a hexagonal shape with the Cu (111) orientation The growth of Cu on ZnO can also be predicted using COMB by calculating the free energ ies of Cu and ZnO surfaces and the interfacial energy of Cu/ZnO, which are difficult to measure experimentally In particular, t hree dimensional film growth is expected when and the two dimensional film growth is expected when where are the free surface energies of Cu and ZnO, respectively, and is the interfacial energy of Cu/ZnO. In addition, based on the charge analysis predicted by COMB, Cu atoms are predicted to have a positive partial charge ( q avg = +0.28 e ) when they were deposited on ZnO but a near neutral charge ( q avg = +0.07 e ) when deposited on the Cu surface Thi s finding indicates that the properties of the initial deposited Cu atoms are highly dependent on the nature of the surface 3 5 Summary T he classical MD simulation results reported here make use of a new, third generation COMB potential for ZnO. They predict that two dimensional (2D) clusters of Cu form on ZnO following the initial deposition of Cu 6 with incident energies of 1.25 eV/atom unt il the coverage of Cu is above 0.4 ML. With increasing coverage, the
63 2D Cu clusters spread to cover more of the clean ZnO surface. However three dimensional (3D) cluster formation is also predicted. In particular, beyond a certain coverage, the transformat ion of 2D clusters into 3D clusters occurs more quickly than the spreading of 2D clusters to fill in the gaps between them. This implies that 3D growth is dominant at higher Cu coverage. The results thus elucidate the ways in which Cu growth transitions fr om layer by layer to three dimensional as the coverage increases. The simulations also consider t he kinetic effect s on Cu clusters and thin film s growth in terms of incident energy and ZnO surface temperature The simulations predict that incident energy i s a critical factor for controlling the final structure of the thin film, while Cu Zn surface mixing may result at high incident energies On the other hand, although increasing the surface temperature does not show a significant influence on surface morph ology as the effect of incident energy did, small and random 3D clusters with compact 2D clusters are found when the surface temperature is maintained a t 800 K. Finally, we compare Cu growth on different supports of Cu (111) and ZnO The simulations clearly indicate that deposited Cu clusters form laminar thin films on Cu (111) and clusters on ZnO In addition, the charge on the Cu atoms is highly dependent on the nature of support ing surfaces. Specificall y, Cu atoms have partial positive charges when they are deposited on ZnO but near neutral charges when deposited on the Cu surface These results are expected to guide the interpretation of experimental data and optimization of supported metal systems for future catalyst design.
64 Table 3 1 Comparison of the properties of ZnO wurtzite phase predicted by the COMB2 and COMB3 potentials with values obtained from experiments and DFT calculations. Properties Exp. DFT COMB2 COMB3 a 0 () 3.242 3.292 3.249 3.267 c 0 () 5.187 5.293 5.213 5.189 c/a 1.600 1.608 1.604 1.582 E c (eV/ZnO) 7.520 7.692 7.440 7.524 B (GPa) 136 134 147 138 C 11 (GPa) 207 230 252 206 C 12 (GPa) 121 82 129 114 C 13 (GPa) 106 64 98 97 C 33 (GPa) 210 247 267 217 C 44 (GPa) 43 75 37 44 C 66 (GPa) 74 57 47 (J/m 2 ) 2.30 2.12 2.21 (J/m 2 ) 2.50 2.31 2.33 (J/m 2 ) 4.30 2.02 2.61 E (ZnS Wurtzite)(eV/ZnO) 0.013 0.008 0.021 E (NaCl Wurtzite)(eV/ZnO) 0.237 0.217 0.903 E (CsCl Wurtzite)(eV/ZnO) 1.358 0.958 1.225 Table 3 2 Properties of Cu Zn alloys predicted by COMB3 compared to those from experiments and DFT calculations. The reference states for the enthalpies of formation and metallic fcc Cu and hcp Zn, which are 3.725 and 1.272 eV/atom, respectively. + Our DFT calculations Phases Properties Exp. DFT + COMB3 Cu 5 Zn 8 a 0 () 8.86 105 8.88 8.88 f (kJ/mol) 126 138 B (GPa) 123 105 156 CuZn a 0 () 3.34 3.36 f (kJ/mol) 21 48 Cu 3 Zn a 0 () 2.97 3.01 f (kJ/mol) 16 36
65 Table 3 3. The adsorption energy (eV) of a single Cu adatom at the three different adsorption sites on the ZnO surface illustrated in Figure 2 as predicted by COMB3 and DFT 107 E ads (eV) Site 1 Site 2 Site 3 DFT 1.52 (0.32 e ) 1.96 (0.55 e ) 1.65 (0.40 e ) COMB3 1.30 (0.55 e ) 1.76 (0.63 e ) 1.47 (0.61 e ) Table 3 4. Predicted migration barriers of a Cu adatom on the ZnO surface found using the dimer method with COMB3 and compared with DFT results 107 E barrier (eV) DFT 0.40 0.30 COMB3 0.76 0.68
66 Figure 3 1 Top view of three different adsorption sites of Cu atoms on the ZnO surface. The Zn (grey) and O (red) are distinguished by color and the atoms in the first (large) and the second (small) layers are distinguished by size. The numbered yellow spheres represent the possible adsorption sites for the Cu adatoms. Figure 3 2 The binding energy (eV) of Cu n (n = 2 to 5) clusters on ZnO predicted by DFT calculations and COMB poten tial. The trendline and R squared value (0.983) are given in the figure.
67 A B Figure 3 3 Relaxed structure of Cu 5 clusters supported by the ZnO surface A) Flat. B) Stack. (At left) Y ellow atoms are Cu, red atoms are O, and grey atoms are Zn). (At right) Atoms are color coded by the charges with the color indicating the charge values.
68 Figure 3 4 Simulation model for deposition of Cu clusters on the ZnO surface
69 A B Figure 3 5 MD simulation o f the morphology of Cu clusters following the deposition of 1ML of Cu atoms. A) T he Cu, O, and Zn atoms are represented by yellow, red, and grey spheres, respectively B) Top down snapshot of just the Cu atoms, which are color coded by height according to the scale shown. Figure 3 6 Snapshots of just the Cu atoms at the indicated coverage (0.2 1ML) following deposition on ZnO The Cu atoms are color coded by height according to the scale shown.
70 A B Figure 3 7 The deposition processes of Cu clusters landing on ZnO surface. A) The deposited Cu cluster lands on the edge of an existing 2D cluster and Cu atoms s tep down to the lower terrace. B) The deposited Cu cluster lands around the center of an existing 2D cluster and starts to form the second layer. Figure 3 8 Schematic of Ehrlich Schwoebel barrier ( ES ) for an adatom on the step edge to descend to the lower terrace in addition to the diffusion barrier ( E diff ).
71 A B C Figure 3 9. Top d own snapshots from MD simula tion of the surface morphology of 1ML of Cu deposited on ZnO w ith different incident energies. A) 0.30, B) 1.25, and C) 2.80 eV/atom.
72 A B C Figure 3 1 0 Top down snapshots from MD simulation of the morphology of 1 ML Cu deposited on ZnO maintained at the indicated temperatures. A) 300 K, B) 600 K, C) 800 K. The Cu atoms are color coded by height according to the scale shown.
73 A B C D Figure 3 11. Snapotshot from an MD simulation of the morphology of Cu clusters on ZnO at the indicated coverage A) 0.5 ML, B) 1.0 ML, C) 1.5 ML, D) 2.0 ML. The Cu atoms are color coded by height according to the scale shown.
74 A B C D Figure 3 12. Snapotshot from an MD simulation of the morphology of Cu clusters on Cu(111) at the indicated coverage. A) 0.5 ML, B) 1.0 ML, C) 1.5 ML, D) 2.0 ML. The Cu atoms are color coded by hei ght according to the scale shown.
75 A B Figure 3 13 Top down snapshots from the MD simulations of 2ML Cu deposited on Cu(111) and ZnO A) Cu(111), B) ZnO The Cu atoms are color coded by height according to the scale shown.
76 CHAPTER 4 DEVELOPMENT OF COMB3 TI/TI O 2 POTENTIAL AND ITS APPLICATONS ON CU CLUSTERS ON TIO 2 (110) SURFACES Titania (TiO 2 ) has attracted much attention due to its diverse applications in pigments, gas sensors, solar cells, photo catalysts and heterogeneous catalysts. 118 It is well known to crystallize into three n atural polymorphs: rutile, anatase and brookite. In addition to these naturally occurring phases, several high pressure phases of titania have been found experimentally or suggested by theoretical calculations, 119 including PbO 2 ), cotunnite (PbCl 2 ), baddeleyite (ZrO 2 ), fluorite (CaF 2 ), pyrite (FeS 2 ), hollandite (MnO 2 ), and ramsdellite (MnO 2 ). Several empirical, fixed charge rigid ion potentials have been developed for atomistic and molecular dynamics simulations of TiO 2 systems 120 121 The form al charge model by Collins et al. 120 and the partial charge model by Matsui et al 122 describe the relative phase stabilities of titania polymorphs and low Miller index rutile surfaces quite well. 119 Following the theories of Sanderson, 123 Parr et al 124 and Mortier et al 125 for charge equilibration ( QEq), Rapp and Goddard introduced the screening the Coulomb interactions that turned the Hessian positive definite for most systems at regular densities. 126 S ome of empirical potentials adopting this QEq scheme have also been developed for TiO 2 including the Morse stretch ( MS Q) 119 potential, the second moment Buckingham (SMB Q) 127 potential and the secon d moment tight binding (SMTB Q) 128 potential. Most recently, a ReaxFF Ti/TiO 2 potential has been developed for describing chemical reactions in Ti/O/Cl/H systems 129 130 In this work, the third generation COMB potential is param eterized for titanium metal and titania and used in classical molecular dynamics (MD) simulations in combination with a previously parameterized potential for Cu to investigate the
77 properties of copper clusters on TiO 2 (110). The rest of the secyion is orga nized as follows. In section 4 1 we introduce the COMB potential functions for the Ti and TiO 2 system s including the fitting procedure properties of the bulk polymorphs and relative stability of Ti and TiO 2 rutile surfaces In s ection 4 2 we demonstrate t he transferability of the COMB potential by applying it to the Cu/TiO 2 system In particular, the adsorption behavior of Cu clusters of different sizes on perfect, reduced and oxidized TiO 2 (110) surfaces are investigated in classical MD simul ations. The findings are compared to the results of density functional theory (DFT) calculations and published experimental results. Lastly, in s ection 4 3 we summarize the major conclusions of this work. 4. 1 Potential Development and Properties of Ti and TiO 2 4.1.1 Fitting Procedure We apply an algorithm that fits structural and energetic information for multiple phases for the system of interest. The procedure determines the potential parameter set using a least squares method to minimize a penalty funct ion, f(p) for each trial parameter set p The penalty function is determined for all included phases by calculating the variation of lattice parameters, elastic constants, and phase stability with corresponding weighting factors. Here, the hexagonal close d packed ( hcp ) and rutile structures are the main focus phase for the Ti and TiO 2 systems, respectively. Therefore, the largest weighting factors are assigned to the properties of these two phases. Wherever possible we use first principles data from the literature, and when this is not available we perform our own DFT calculations to complete the data set for fitting and comparison. Our first principles calculations are carried out with the Vienna ab initio Program (VASP) 59 62 program using the gene ralized gradient approximation (GGA) 131
78 and the Perdew Burke Ernzerhof ( PBE) 43 exchange correlation functional. The calculations utilize plane wave basis sets with a 500 eV energy cutoff, a 444 Monkhorst Pack k point mesh 132 and projector augmented wave (PAW) pseudopotentials 133 134 for Ti ([Ne] core with 3s 2 3p 6 4s 2 3d 2 electrons) and O ([He] core with 2s 2 2p 4 electrons). The convergence criteria are set at 1.0 10 5 eV and 1.0 10 3 1 for energies and forces, respectively. The procedure used for parameter fitting follows that outlined by Liang et al. 135 In particular, element type parameters are fit for pure elements (Ti and O) and bond type parameters are fit for TiO 2 systems. The fitting structures for metallic Ti include the hcp face centered cubic ( fcc ), body centered cubic ( bcc ), simple cubic, and diamond phases. The element type parameters of oxygen, which are kept the same throughout the third generation of COMB to ensure full compatibility and transferability of the different potentials, are fit to dimer and trimer oxyge n molecules. For TiO 2 binary systems, we include 10 different AB 2 phases in the fitting training set, following the work of Swamy et al 119 These include rutile (space group: P4 2 /mnm ), an atase ( I4 1 / amd ), and brookite ( Pbca ) in addition to the high pressure phases columbite ( Pbcn ), cotunnite ( Pnma ), baddeleyite ( P2 1 /c ), and fluorite ( Fm 3m ) and the hypothetic phases of pyrite ( Pa 3m ), hollandite ( I4 / m ), and ramsdellite ( Pbnm ). The atomic a nd electrostatic parameters for Ti and O while the bond dependent parameters for Ti O and O Ti are provided in A 1,A 2, and A 4 4.1.2 Predicted Ti Properties The static properties of metallic Ti are reported in this section, where the structures are geometrically optimized until the deviation of total energy and force on
79 each atom are smaller than 1x10 5 eV and 1x10 2 eV/, respectively. Table 4 1 compares the p hysical properties of the hcp phase of Ti predicted by COMB with values from published experimental data, DFT calculations, and a modified embedded atom method (MEAM) potential. The properties marked with an asterisk denote fitted properties, while all oth ers are predicted properties. The COMB potential successfully reproduces the phase stability of Ti, with hcp the ground state phase, and fcc and bcc having higher energies. The fitted energy differences between fcc hcp and bcc hcp are 65 meV/atom and 100 meV/atom, respectively, which are in a good agreement with the DFT values. The lattice parameters of hcp Ti metal fitted by COMB deviate from the experimental values by 0.5% and + 1.1% for a 0 and c 0 respectively. These deviations result in a larger valu e of c/a ratio=1.61 which is larger than the experimental value of 1.587 Most importantly, both are significantly lower than the ideal value of c/a=1.633; it can thus be expected to reproduce the experimentally observed deformation mechanisms. 136 137 T he cohesive energy and elastic constants fitted by COMB reasonably well reproduce the experimental and DFT values. The vacancy and interstitial defect formation energy are qualitatively reproduced by COMB, although the Ti interstitial (octahedral site) formation energy is overestimated. The res ult indicates that t he formation energy for a Ti vacancy is 2 eV smaller than the energy to form a Ti interstitial. Finally, the predictions for surface energies of the three low Miller index surfaces are given in Table III. The COMB potential predicts tha t the (0001) basal plane has a lower surface energy than the surfaces on the and planes. This trend is also predicted by the MEAM
80 potential. However, these predicted trends do not agree with the results of the DFT calculations, which predict the prismatic surface to be the most stable Ti surface. The average predicted surface energies from the COMB potential are ~25% smaller than the average value measured experimentally. 4.1.3 Properties of Rutile TiO 2 and Phase Stability The properties of rutile TiO 2 from the COMB potential, compared to the experimental values, DFT calculations, and three other classical interatomic potentials (Matsui 122 MS Q 119 SMB Q 127 and ReaxFF 129 ) are presented in Table IV. The Matsui and MS Q values are generated with the General Utility Lattice Program (GULP) 138 139 us ing parameters provided in Refs 119 and 122 respectively. For the rutile phase, the COMB potential fits the lattice parameters and cohesive energy in a good agreement with experimental values, where the deviations are all below 1%. Wit h regard to the elastic constants, although the larger overestimation of C 12 leads to a larger bulk modules reproduced by COMB, the reproduced elastic constants of C 11 C 33 C 13 C 33 C 44 and C 66 are within 20% of the experimental and DFT values. One thin g to note is that the violation of the Cauchy relation is observed experimentally. 140 141 The Cauchy relation states that C ijkl =C ikjl so that, for example, C 12 is equal to C 44 for cubic systems and C 12 is equ al to C 66 for tetragonal systems such as rutile. To faithfully reproduce the elastic moduli of TiO 2 COMB should capture this violation so that C 12 C 66 for TiO 2 which as indicated in Table 4 2 it does reasonably well. Table 4 3 lists the energy differences between the rutile and other AB 2 phases predicted by COMB, DFT and other empirical potentials. Both the COMB and Matsui
81 potentials predict the rutile phase to be the most stable, followed by columbite, brookite and anatase. Th e MS Q and SMB Q potentials also predict rutile to be the ground state phase, but followed by brookite and anatase, respectively. Although the COMB potential does not predict exactly the same phase order for all other AB 2 phases relative to DFT, rutile is correctly predicted to be more stable than the high pressure or hypothetical AB 2 phases. This structural stability is important when the properties of rutile phases are determined at high temperatures which may lead to phase transformations. 4.1.4 TiO 2 Low Miller Index Surfaces Single crystalline TiO 2 surfaces are important in a wide range of applications. 118 An important aspect of assessing the quality of a TiO 2 empirical potential is its ability to describe surface properties. We therefore investigate the capability of the developed COMB TiO 2 potential by calculating the surface energies of pristine TiO 2 rutile low Miller Index surfaces. Since experimental numbers for these surface energies are difficult to obtain, we largely compare values predicted from COMB to the results of DFT and other empirical potentials. Table 4 4 compares the energ ies of the (110), (100), and (001) surfaces of rutile predicted by the COMB potential, DFT 142 143 and other empirical potentials 119 122 127 144 145 The relaxed (unrelaxed) surface energies predicted by COMB for (110), (100), and (001) are 0.99 (1.09), 1.25 (1.76), and 1.49 (1.87) J/m 2 respectively. The surface energies co mputed with first principles techniques over the last fifteen years are reported by Hallil et al 145 Because different approaches (linear combinations of atomic orbitals, plane waves), software packages, and number of surface slab layers are used in the calculations, the (110) surface energies range from 0.40 to 1.10 J/m 2 Nevertheless, this low Miller index orientation is predicted to be the most stable surface
82 by all methods considered, followed by the (100) and (001) orientations for the TiO 2 rutile phase. Variable charged empirical potentials such as COMB, SMB Q 127 and SMTB Q 145 predict the correct order of surface stability, where SMTB Q does a better job in predicting the value of surface energies relative to DFT. Kim et al. 129 report ed the surface energies for five different surfaces (t hree for the anatase phase and two for rutile) using ReaxFF, which predicts that (110) is more stable than (101) surfaces in the case of rutile. In contrast to predictions from the above potentials, the MS Q potential predicts an incorrect order for the (1 10) and (100) surface energies. In addition to the surface energy, it is important that the property of surface atoms, such as charge charger, is properly captured as they influence processes such as catalysis and thin film growth. The charge transfer defined here is the charge variation between the atomic charge of a surface atom and its charge in the bulk. The reference bulk charges for Ti and O in COMB are +1.905 and 0.953, respectively. The surface atoms selected for charge transfer determination on (110), (100), and (001) surfaces a re labeled in Figure 4 1. Table 4 5 presents the charge transfer determined by the COMB potential compared with the published results of Mulliken charges from DFT calculations with the CRYSTAL06 code. 145 For the (110) and (100) surfaces, the predictions of COMB are in a good agreement with the results of DFT, particul arly for the sign of charge transfer. Both show only bridge oxygen atoms O (1) becoming more positively charged with respect to the bulk, although the charge transfer on O (1) of TiO 2 (110) is overestimated by COMB. In the case of the (001) surface, charge tr ansfer also takes place primarily on
83 the most undercoordinated surface atoms, Ti (1) and O (1) In particular, COMB predicts that the four fold coordinated Ti (1) atom become more negatively charged with respect to the bulk (six fold) and other five fold Ti s urface atoms of the (110) and (100) surfaces, which is consistent with the results of DFT calculations. This result demonstrates that COMB is capable of describing correctly the charge of atoms under different environments, which is an important criterion for the modeling heterogeneous catalysis. 4.2 Heterogeneous Interfacial System: Cu n /TiO 2 (110) 4.2.1 Cu Clusters on Stoichiometric TiO 2 (110) Understanding the growth mode of metals on oxide surfaces is important for the design of electronic devices, catalysts, and sensors. 1 46 The main strength of the COMB potential is its transferability that allows the combination of established potentials for modeling heterogeneous systems, such as metal/metal oxide systems. In this section, we consider the energetics of Cu clusters sup ported by rutile TiO 2 (110), which is of particular interest for heterogeneous catalysis 147 To characterize the relative stability of diff erent sized copper clusters, the binding energy ( E b ) per atom is estimated using the following standard equation: (4 1) where n is the number of Cu atoms, E ( Cu n / Ti O 2 ) is the total energy of Cu n adsorbed on the surface, and E ( Cu n ) and E ( Ti O 2 ) are the energies of the Cu n clusters and the optimized TiO 2 surface slab, respectively. The preferred adsorption site and energetics for single Cu atoms on this surface was inve stigated first. A Cu atom was placed above the bridging O atom (labeled O (1) in Fig. 4 1), the five fold coordinated Ti atom (labeled O (1) in Fig. 4 1), and in between two
84 bridging O atoms. The COMB potential predicts the Cu binding energies on these sites to be 2.32, 0.15, and 3.25 eV, respectively. Thus, COMB predicts that the most favorable adsorption site for Cu is in between two bridging O atoms, which is also predicted by Giordano et al 148 and Pillay et al 149 using DFT calculations and observed by experiments. 150 This agreement gives confidence regarding the interaction of Cu cl usters with this TiO 2 surface, something that is more computationally intensive to investigate using first principles methods. The configurations of the Cu clusters considered are illustrated in Fig ure 4 2. Cu 13 Cu 55 and Cu 147 are constructed with a cub o octaherdron shape, while Cu 38 is constructed with a truncated octahedron shape. Based on the binding energy, Cu 147 interacts most strongly with the surface than the other clusters considered. This higher binding energy is thought to result from the forma tion of more Cu O bonds across the Cu cluster/TiO 2 interface. The charge distributions of the different sized Cu clusters are illustrated in Fig ure 4 3. Charge transfer is predicted to occur between the clusters and TiO 2 surface by the COMB potential that results in positively charged Cu clusters in agreement with DFT calculations. 148 The average charge of Cu 13 Cu 55 and Cu 147 clusters predicted by COMB are 0.33, 0.21, 0.20, a nd 0.15 respectively. 4.2.2 Cu 147 Clusters on Non Stoichiometric TiO 2 (110) Previous experimental and theoretical studies have reported that the growth modes and optimal size of metal clusters deposited on TiO 2 depends on the details of the surface structure. 151 154 To quantify the influence of surface defects on the binding energy of Cu clusters, a reduced rutile (110) surface is built by removing bridging o xygen surface atoms from a perfectly stoichiometric surface slab. Similarly, an oxidized
85 TiO 2 (110) surface is created by adding oxygen atoms on top of five fold surface titanium atoms. The largest binding energy predicted for the systems considered here i s calculated for the Cu 147 cluster on oxidized TiO 2 (110). This enhanced bonding is in agreement with both experimentally observed findings and DFT calculations for clusters of Au 151 153 and Ag 152 153 on the same oxidized surface. In addition, through charge distribution analysis of Cu 55 on non stoichiometric TiO 2 (110), shown in Fig ure 4 4, it is clear that while the charge distrib ution of external Cu cluster atoms are similar as the oxidation state of the surface changes, the number of positively charged Cu atoms at interfaces increases when the clusters are supported on the oxidized surface. The positively charged Cu atom correspo nds to a metal oxide bond between Cu and O that is responsible for the higher bonding energy to the surface. The agreement between the predictions of the COMB potential, experimental data, and DFT calculations indicates that third generation COMB potentials for Ti/TiO 2 and Cu/Cu 2 O may be used to model Cu/TiO 2 interfaces. However, to capture more complex phenomena, such as Cu cluster or thin film growth, additional modification of the potential fitting database to emphasize quantities such as diffus ion barriers may be necessary. 4.3 Summary The work described here summarizes the development of an empirical, variable charge potential for Ti and TiO 2 systems within the third generation COMB potential framework. We also carried out first principles DF T calculations to determine the phase order of TiO 2 polymorphs. The developed potential is fitted to, and captures most of, the structural, energetic and mechanical properties of Ti metal and TiO 2 rutile phases. In addition, rutile is correctly predicted t o be more stable than the high pressure or
86 hypothetical AB 2 phases. The potential generally reproduces the surface stability of three low Miller index rutile surfaces although their energies are overestimated. Combining the Ti/TiO 2 potential developed in this work with the previously developed COMB potential for the Cu/Cu 2 O 155 system, we d emonstrate that the adsorption energies of Cu atoms and clusters on the TiO 2 (110) surface strongly depend on Cu O interactions. Higher binding energies result from a larger number of Cu O bonds formed at the Cu/TiO 2 interface. Therefore, in this work, COM B predicts that Cu 147 has the largest binding energy compared to the other clusters considered. Charge transfer is predicted to take place between Cu clusters and the TiO 2 surface that results in positively charged clusters. Additionally, COMB predicts an enhanced bonding between Cu clusters and the oxidized surface, which is consistent with findings from both experimental observations and DFT calculations for late transition metals (Au and Ag) on oxidized TiO 2 (110). Thus, the third generation COMB potent ial for Ti/TiO 2 is shown to successfully capture many of the important trends for these materials and be transferrable when combined with previously developed potentials to model heterogeneous interfacial systems.
87 Table 4 1. Properties of Ti reproduced or predicted by the COMB3 potential developed in this work, and compared to experimental results, DFT calculations and the MEAM potential. Properties Exp. DFT + COMB MEAM + a 0 () 2.951 2.948 2.935 2.931 c 0 () 4.684 4.667 4.705 4.678 c 0 /a 0 1.587 1.583 1.602 1.596 E c (eV/atom) 4.844 5.171 4.843 4.831 C 11 (GPa) 176 172 156 174 C 12 (GPa) 87 82 86 95 C 13 (GPa) 68 75 80 72 C 33 (GPa) 191 190 170 190 C 44 (GPa) 51 45 31 58 C 66 (GPa) 28 ++ 36 B (GPa) 110 111 104 113 G (GPa) 52 47 36 55 E (fcc hcp) (meV/atom) 60 58 65 39 E (bcc hcp) (meV/atom) 70 108 100 111 Vacancy (eV/defect) 2.03 2.34 2.24 Interstitial (eV/defect) 2.58 4.39 2.64 (mJ/m 2 ) 1920 1939 ++ 1341 1474 (mJ/m 2 ) 1920 2451 ++ 1405 1554 (mJ/m 2 ) 1920 1875 ++ 1677 1682 Denotes fitted properties for the COMB potential + DFT and MEAM values from Ref.[ 156 158 ] unless otherwise specified ++ Our DFT calculation
88 Table 4 2. Properties of the rutile, anatase, and brookite phases calculated by the COMB potential. Th e properties are compared to experiments, DFT calculations (values without references are our calculations), and three other empirical potentials. Rutile TiO 2 Properties Exp. DFT COMB SMB Q MS Q Matsui h a 0 () 4.594 4.646 4.562 4.581 4.587 4.493 c 0 () 2.959 2.970 2.967 2.966 2.958 3.008 E c (eV/TiO 2 ) 19.90 26.95 19.19 19.90 16.49 39.80 B (GPa) 218 a 211 b 266 228 229 237 G (GPa) 124 a 136 b 125 87 116 C 11 (GPa) 268 a 270 b 318 293 294 322 C 12 (GPa) 175 a 172 b 257 203 202 230 C 13 (GPa) 147 a 147 b 182 164 168 147 C 33 (GPa) 484 a 468 b 516 400 423 444 C 44 (GPa) 124 a 116 b 123 128 96 123 C 66 (GPa) 190 a 216 b 204 183 190 226 q Ti ( e ) 2.26 1.91 2.51 1.15 2.20 Anatase TiO 2 a 0 () 3.785 3.806 3.849 3.825 3.850 3.770 c 0 () 9.512 9.714 9.135 9.030 9.060 9.568 E c (eV/TiO 2 ) 27.03 19.06 19.84 8.19 39.49 B (GPa) 59 178 c 180 234 220 176 176 G (GPa) 57 126 66 59 q Ti ( e ) 2.24 1.88 1.12 2.20 Brookite TiO 2 a 0 () 9.184 d 9.368 9.259 9.113 9.146 b 0 () 5.449 d 5.424 5.444 5.449 5.389 c 0 () 5.145 d 5.046 5.229 5.170 5.144 E c (eV/TiO 2 ) 19.11 19.65 8.21 39.62 B (GPa) 255 e 261 227 211 200 G (GPa) 124 83 92 q Ti ( e ) 1.89 1.14 2.20 Denotes fitted properties for COMB potential a Reference 159 h Reference b Reference 160 g Reference c Reference 161 d Reference 162 e Reference 163 f Reference
89 Table 4 3 Cohesive energies of rutile titania (eV/TiO 2 unit) and that relative to rutile titania for other polymorphs (AB 2 phases) calculated from COMB potential for TiO 2 developed in this work in comparison with DFT calculations and other empirical potentials. Polymorph (space group) Exp. DFT COMB3 MS Q b Matsui c Normal phases Rutile (P4 2 /mnm) 26.949 19.189 8.248 39.798 19.90 Anatase (I4 1 /amd) 0.035 a 0.083 0.129 0.057 0.304 0.060 Brookite (Pbca) 0.008 a 0.046 0.080 0.037 0.178 0.250 High pressure phases Columbite (Pbcn) 0.006 0.010 0.074 0.035 Cotunnite (Pnma) 0.107 0.739 0.553 0.914 Baddeleyite (P2 1 /c) 0.066 0.667 0.515 Fluorite (Fm 3m) 0.775 1.118 1.020 0.580 Hypothetical phases Pyrite (Fm 3m) 0.253 Hollandite (I4/m) 0.275 Ramsdellite (Pnma) 0.263 a Reference 164 b Reference c Reference d Reference
90 Table 4 4. Low index s urface energies (J/m 2 ) of rutile TiO 2 from the ab initio approache s and empirical potentials. Method Surface formation energy (J/m 2 ) (110) (100) (001) ab initio PBE (Crystal03) 143 0.42 0.69 1.39 LDA (Crystal03) 143 0.90 1.20 1.88 B3LYP (Crystal03) 142 1.10 Empirical potentials Fix charged CFR 119 14 4 0.95 1.24 2.20 MA 119 122 1.77 2.07 2.40 Variable charged MS Q 2 1.40 1.32 1.85 SMB Q 127 0.60 0.85 1.74 SMTB Q 145 0.42 0.49 1.26 COMB 0.99 1.25 1. 4 9
91 Table 4 5. Charge variation (Q Q bulk ) (in e ) of surface atoms (labeled in Figure 4 1). The charges of Ti and O in bulk rutile are1.905 and 0.953, respectively Q Q bulk CN DFT COMB TiO 2 (110) O(1) 2 0.18 0. 28 O(2) 3 0.04 0.0 8 Ti(1) 5 0.03 0.02 Ti(2) 6 0.06 0.08 TiO 2 (100) O(1) 2 0.1 8 0.2 0 O (2) 3 0.02 0. 09 O(3) 3 0.0 1 0.0 4 Ti(1 ) 5 0. 09 0.0 5 TiO 2 (001) Ti(1) 4 0.15 0.09 Ti(2) 6 0.05 0.02 O(1) 2 0.05 0.11 O(2) 3 0.04 0.05
92 A B C Figure 4 1 Model of the rutile TiO 2 surface s A) (110), B) (100), and C) (001) surfaces. Ti and O atoms are represented by big (grey) and small (red) spheres. The surface atoms for charge transfer calculations are labeled by element type and Arabic numbers.
93 Figure 4 2 Optimized geometries of Cu 13 Cu 38 Cu 55 and Cu 147 clusters. Figure 4 3 Charge distribution of relaxed Cu clusters on TiO 2 (110) A ) Cu 13 B ) Cu 38 C ) Cu 55 and D ) Cu 147 Large atoms are Ti, medium sized atoms are Cu and small atoms are oxygen.
94 Figure 4 4 Charge distribution of relaxed Cu 55 clusters on the stoichiometric and defected TiO 2 (110) surfaces. A) R educed, B) stoichiometric, and C) oxidized TiO 2 (110) surfaces Large atoms are Ti, medium sized atoms are Cu and small atoms are oxygen.
95 CHAPTER 5 DEVELOPMENT OF COMB3 TI/ T IN POTENTIAL AND ITS APPLICATONS ON TI/TIN INTERAFACES AND THE DEPOSTION OF ATOMIC OXYGEN ON TIN(001) SURFACE Titanium nitride (TiN) thin film s are notable for their distinctive properties, which include extreme hardness (2000 kg / mm 2 ), high melting temperat ure (2950 C), and a small thermal expansion coefficient (9.4 x 10 6 C 1 ). 38 The y have consequently been widely used as wear resistant coatings on cutting and grinding tools 3 165 166 corrosion resistant coverings on mechanical parts 167 168 and diffusion barriers in microelectronic devices 169 More recently, TiN thin films have been used in medical implants 170 171 because of their biocompatibility with bony tissue. All of these applications are highl y dependent on the nature of the interfac e structure between TiN and the underlying substrate. Therefore, achieving a better understanding of the properties of these highly heterogeneous interfaces across length scales will allow us to better enhance the u tilization of TiN thin films. In the case of the TiN system, the bulk and surface properties of TiN have been studied by Marlo et al 172 using density functional theory (DFT) with various exchange correlation functions However, there are only a few empirical potentials for TIN in the literature. The first interatomic potential description of TiN was reported by Hernndez et al 173 in the form of a rigid ion model. However, no information was provided on how this potential can reproduce the fundamental physical properties of TiN system. Later, an exten ded Tersoff type potential of TiN was developed by Iwasaki et al 174 In this case, although the details of the potential formalism and parameters are provided, the structural and elastic properties of TiN were not discussed Kim et al 175 developed a potential ba sed on the second nearest neighbor modified embedded atom method
96 (2NN MEAM). They provid e d the parameters and to extensive ly discuss ed its ability to predict the properties of TiN. Later, the parameters of Kim et al w ere modified slightly by Yu et al 156 for a better description of physical properties of Ti 2 N. Most recently, Xu et al 157 combined the 2NN MEAM potential with Ziegle r Biersack Littmark (ZBL) repulsive potential to examine the growth of TiN films on the TiN(001) surface in MD simulation; the results indicated that the TiN film grows via layer by layer growth mode. The purpose of the present work is to develop an interatomic, third generation COMB potential for the Ti N system compatible with the developing potential of Ti O and N O systems This thus enables the simulation of complex Ti N O phenomena, such as oxidation of TiN. 5 .1 Potential development and Propert ies of Ti, N 2 and TiN 5 .1.1 Fitting Procedure The element type parameters are fit for pure elements (Ti and N) and bond type parameters for binary systems (TiN). In the present work, the element type parameters of Ti are implanted from the established COM B potential for Ti/TiO 2 system, 176 where the fitting structures for metallic Ti include the hexagonal c losed packed ( hcp ), face centered cubic ( fcc ), body centered cubic ( bcc ), simple cubic, and diamond phases. For N, the binding energy and bond length of N 2 molecule are fitted to experimental values. For the bond type parameters, to ensure the phase stability of the ground state, the fitting structures for TiN system include rocksalt ground state ( B1 space group: ), wurtzite ( B4 space group: ), zincblende ( B3 space group: ), and cesium chloride ( B2 space group: ). Because Ti 2 N c an form stable intermediate phases according to the Ti N phase diagram, we also consider four fitting structures for
97 the Ti 2 N system, which are rutile (space group: ), which is the ground state structure of Ti 2 N, fluorite (space group: ), F 2 N type (space group: ) and Cu 2 O type (space group: ) structures. All the fitted structures considered are illustrated in Figure 5 1. T he physical properties of the TiN rocksalt phase are of the most interest. We therefore include the lattice constants, elastic properties, defect formation energies, and surface energies of this phase in the fitting from published experimental data and DFT calculations. In those cases where reliable data could not be obtained from the literature, we performed our own DFT calculations to complete the data set for fitting and comparison. The DFT calculations were carried out with the Vienna ab initio Simulation Program (VASP) 61 62 177 178 within the generalized gradient approximation (GGA), using the Perdew Burke Ernzerhof (PBE) 43 e xchange correlation functional. Ti3 p 6 3 d 3 4 s 1 and N2 s 2 2 p 3 were taken as the valence electrons A plane wave basis with a 400 eV energy cutoff and a 666 Monkhorst Pack k point mesh 179 180 were used. The convergence criteria for energies and forces we re set at 1x10 5 eV and 1x10 2 1 respectively. For N 2 molecule, the electronic structure calculations were performed using the Gaussian09 181 computational chemistry software package. The geometric optimization was carried out using B3LYP/6 311 G*, 63 and the ionization energies are calculated with the CCSD(t)/cc pVTZ 182 183 level of theory. 5 .1.2 Fitte d Properties of Ti metal and N 2 molecule The static properties of lattice constants for Ti metal are provided in Table 4 1 where the hcp structure is the ground state phase for Ti metal. In the case of the N 2 molecule, the bond length and bonding energy match well the experimental values, 184
98 which are 1.098 for COMB vs. 1.099 in experiment and 9.76eV vs 9.80 eV per molecule. 5 .1.3 Fitt ed and Predicted TiN Phase Stability and Properties Table 5 1 compares the physical properties of TiN polymorphs fitted or predicted by th e COMB potential with values from published experimental data, DFT calculations, and second nearest neighbor modified embedded atom method (2NN MEAM) potential 184 In order to determine the phase stability in TiN and Ti 2 N systems, the enthalpy of formation, H f is fitted. The definition of the enthalpy of formation of a compound is the change of enthalpy from the formation of one mole compound from its constitutive elements in their standard states: hcp Ti metal and an isolated N 2 molecule re spectively. As shown in Table 5 1 t he rocksalt and rutile phase s are correctly fitted as the most stable for TiN and Ti 2 N systems, respectively For the TiN system, which is the phase of primary interest, the COMB potential not only successfully fits the lattice constant and enthalpy of formation for the rocksalt structure to experimental data, but also reproduces the correct phase order of rocksalt, wurtzite, zincblende, and cesium chloride. The fitted energy differences between the rocksalt phase and hypothetical phases of wurtzite, zincblende, and cesium chloride are also in good agreement with the results of DFT calculations. Wi th regard to the elastic properties, COMB does a better job of fitting the values of C 11 C 22 and C 44 than the 2NN MEAM potential, as illustrated in Table 5 1 In particular, the largest deviation of elastic constants with respect to the experimental valu es is 4% for COMB and 12% for 2NN MEAM due to its larger overestimation of C 44
99 In the case of the Ti 2 N system, the rutile phase is correctly fitted to be the ground state structure. However, the agreement between fitted and experimental values for the la ttice constants and the enthalpy of formation for the rutile structured Ti 2 N are only fair. Nevertheless, COMB successfully reproduces that the formation TiN is favored over the formation of Ti 2 N based on the standard state conditions used in the enthalpy of formation calculations. This is important for the atomic scale study of TiN alloys. 5 .1.4 Predicted Defect Formation Energies and Surface Energies of rocksalt TiN The defect formation energies for the Ti vacancy ( V Ti ), N vacancy ( V N ), Ti interstitial ( I Ti ) and N interstitial ( I N ) and the Schottky and Frenkel complexes have been determined for the TiN rocksalt structure The vacancies are generated by removing an atom of Ti or N from the lattice sites, while interstitials are created by placing a Ti or N atom in the tetrahedral position. The Schottky defects are created by the association of a Ti and an N vacancy that are separated by a distance of 4.77 to prevent their direct interaction. Similarly, to prevent the recombination of Frenkel or anti Fre nkel, the vacancy and interstitial are separated by 4.63 Table 5 2 compares the defect formation energy predictions of the COMB potential with published DFT results 173 The defect formation energies E f were calc ulated with the standard e quation ( 5 1) where E def an d E perf ar e the energies of defective and perfect bulk structures, respectively, n i is the number of defects and i is the chemical potential of the atom added or removed. These energies were not part of the fitting dataset of physical properties ; COMB predicts generally larger values of defect formation energies than
100 DFT. Nevertheless, COMB predicts defect formation energies with a fidelity si milar to that of DFT which the nitrogen vacancy and the Schottky defect predicted to be the most stable point defect and defect complex, respectively. This gives confidence that the COMB potential captures much of the correct physics. In addition to the formation of defects, the predictions for the surface stability of TiN were investi gated. To simulate surfaces, a vacuum of 25 is added in the direction normal to the surface within the simulation box, thereby forming a slab with two free surfaces. To facilitate comparison, the surfaces are relaxed until the force on each atom is close d to zero and there is a minimal fluctuation in total energy. For the (001) and (110) surfaces, which are nonpolar, the number of Ti cations and N anions within each layer is identical. Therefore, a stoichiometric slab with two identical surfaces may be cr eated. In contrast, the TiN (111) surface is polar and is constructed with alternating layers of Ti cations and N anions. It therefore is expected that a rocksalt (111) surface will undergo reconstruction 185 to minimize the surface dipole. In this study, a stoichiometric slab with Ti and N terminated surfaces for TiN (111) is created for the surface energy calculation. following standard equation 185 : ( 5 2) where E slab is the total energy of the slab for the surface energy calculation, E bulk is the energy per atom of the perfect bulk structure, n is the total number of atoms in the slab, and A is the total area of surfaces (top and bottom) of the slab. Table 5 3 compares the surface energies predicted b y COMB to the results of DFT calculations and MEAM
101 potentials. The surface energy trend predicted by COMB is (100) < (110) < (111), which is consist ent with the trend predicted by DFT Overall, COMB potential corre ctly reproduces the phase stability, lattice constant, and elastic constant of TiN rocksalt structure compared to the experimental value and DFT calculations. However, the underestimation for the surface energy of the TiN(111) polar surface needs a further improvement. 5 .2 Predicted Dynamical Results The potentials fitting performed using calculations carried out at 0 K To verify that there are no unphysical p hase transformation during MD simulations at high temperatures the rocksalt phase TiN is heated to 3500 K and then cooled to 300 K 186 The simulation cell has dimensions of 10 a 10 a 10 a ( a = 4 242 lattice constant of TiN) and consists of 8000 at oms in the stoichiometric rocksalt structure. All of the simulations are performed using the LAMMPS software. The system is initially equilibrated at 2000 K for 100 ps in the constant temperature and pressure (NPT) ensemble until the total energy variation and total pressure of the simulation box is below the values of convergence which are both at 1x10 4 eV and 1x10 2 1 Next several successive equilibration calculations are carried out by increasing the temperature at intervals of 100 K up to 3 5 00 K At each temperature, the system is equilibrated for 40 ps Fig ure 5 3(a) shows the change of the crystal configuration as a function of temperature. At 1000 and 2000 K, the rocksalt structure of TiN is stable. At 3000 K, it is clear to see that Ti and N atoms have a larger thermal vibration but are still in their original lattice positions. Finally, at 3500 K, all of the atom positions abruptly randomize, resulting in a liquid phase. On the basis of the radial distribution function that is illustrated in Fig ure 5 3(b), the rocksalt structures of
102 TiN are maintained up to 3000 K. As for the prediction of melting point by MD 76 187 because the MD structure is a defect free crystal, it is expected that the prediction will be higher than the experimental values where the samp le always has defects that facilitate the transition from the crystalline to the liquid phase. Therefore, it is only possible to indicate that the COMB potential predicts a m elting point around 3000 3500 K on the basis of these simulations, which is in excellent agreement with the experimental melting point of 3203 K. 188 For the quenching process, t he structure was equilibrated by follow ing the heating procedure back down to 300 K from 3000 K Our simulation result indicates that the heated structure returns to the rocksalt phase. 5 .3 A pplications of COMB 3 TiN P otential 5 .3.1 Adhesion of fcc Ti/TiN( 001) Interfaces Ti/TiN multilayered structures are commonly used as barriers in u ltra large s cale i ntegration (ULSI) devices for Al metallization 3 or coated on steel to provide corrosion resistance 189 and enhance tribolog ical properties 190 co mpared to single layered hard coatings Because the performance in the above applications is highly dependent on the nature of the interfac e between Ti and TiN, an atomic scale approach provides important insights into the nature of the bonding at metal/ce ramic interfacial structures. Here, to demonstrate the ability of COMB3 TiN potential for describing the different bonding environments in a seamless manner the interfacial energy between Ti and TiN was determined. Based on the x ray diffraction analysis, TiN thin films with (001) orientated surface were observed when they were sputter deposited on the glass substrate. 191 Since experimentally Ti thin films are epitaxially grown on TiN and there are only small differences in the lattice constants (~ 3.3%) between fcc Ti metal (4.10 )
103 and rocksalt TiN (4.24 ), a coherent interface of fcc Ti(001) || TiN(001) is expected. Here, three coherent (001) interfaces were considered between fcc Ti and rocksalt TiN, as shown schematically in Figure 5 4 (a). For cases (1) and (2), the Ti metal interact s directly with Ti atoms and N atoms of TiN respectively. For case (3), the Ti atoms in the metal are on the bridge site of TiN and only interact with Ti atoms in TiN. The work of adhesion, W was used to quantify the stability of in terfaces and was calculated from the standard equation : where E Ti and E Ti N are the total energy of relaxed fixed volume Ti and TiN slabs, respective, E Ti /TiN is the total energy of the relaxed interfacial structure, and A is the interface area. All the structures were equilibrated at 300 K, prior to quenching at T=0 K. The calculated works of adhesion for these three interfaces are 4.45 J/m 2 0.93 J/m 2 and 2.83 J/m 2 respectively. The corresponding values from DFT calculations 192 are 3.56 J/m 2 0.43 J/m 2 and 1.66 J/m 2 respectively. The results indicate that COMB potentials qualitatively predict the trends compared to DFT calculations. In particular, both COMB and DFT predict that the interface of Ti metal interacting with N atoms in TiN is stronger than that of Ti metal interacting with Ti atoms in TiN, which can be attributed to the stronger bonding associated with Ti N bonding relative to Ti Ti bonding. In addition, the atomic charge predicted by COMB for these interfaces, shown in Figure 5 4 (b), indicates that the Ti and N atoms in TiN have more positive and negative charge when Ti atoms in Ti on t he top sites of N or bridges site of Ti in TiN.
104 5 .3.2 Atomic Oxygen deposition on TiN(001) surface The adsorption of oxygen atoms on the TiN surface plays a crucial role in the TiN oxidation processes Here, we perform two MD simulation s to examine the adsorption of atomic and molecular oxygen on the TiN (001) surface through deposition. In the case of atomic oxygen deposition, sixty oxygen atoms are deposited with an incident energy of 5 eV normal to the surface in spatially random loc ations. T he dimensions of the TiN surface slab are 42 42 63 where the bottom layer is fixed in their bulk lattice sites to prevent the slab from translating as a result of the deposition process. The middle eight layers of TiN are thermostat ed u sing a Langevin thermostat applied to dissipate the excess energy from the deposited oxygen atoms and to maintain a constant surface temperature. The active region comprises the top six layers of the slab ; these evolve under unconstrained equations of moti on The surface slab was equilibrated at 300 K for 10 ps prior to the deposition Fig ure 5 5 shows the snapshots after thirty and sixty oxygen atoms have been deposited, respectively. A majority 70% of the deposited O atoms stick to the surface. Interesti ngly, w hen the COMB potential is used to calculate the static adsorption energy for atomic oxygen on this surface the large st adsorption energ y is predicted for the a top site on the N atom, which contradicts the results of DFT calculations 193 However the larger cutoff distance of Ti O (3.2 ) compared to N O (2.3 ) within the potential leads the oxygen atoms to preferably adsor b on top of the Ti atoms in the MD simulations, in agreement with the r esults of DFT calculations The simulations further predict that 1 0% of the incident oxygen atoms rebound from the surface and 20% of the incident oxygen atoms form O 2 molecule s in the vacuum above the surface
105 In the case of the simulation of O 2 deposit ion, a single oxygen molecule is given an given incident energy of 5 eV per molecule and placed 5.0 above a TiN (001) surface slab with the same dimensions and characteristics as described above for the atomic O deposition case. The oxygen molecule binds initially to a Ti atom. Subsequently, it moves to a bridge position over two Ti atoms and then dissociates, as illustrated in Fig ure 5 6. These results are consistent with the findings of a study by Piscanec et al. 194 who analyzed the early oxidation stages of TiN(100) by means of first principles molecular dynamics (FPMD). The bridge Ti atoms will be the preferred adsorption site for the oxygen molecule. Thus, through a combination of characteristics these examples demonstrate the fact that the third generation COMB potential may be used to examine Ti N oxidation via MD simulations. 5 .4 Summary In summa ry an empirical, variable charge potential for TiN systems was successfully developed based on the COMB3 potential framework. The potential captures most of the physical properties of TiN rocksalt phase. Because the TiN is the primary interest in the work it is acceptable that the fitted lattice constants and enthalpy of formation for Ti 2 N rutile phase are only fair compared to experimental values. Nevertheless, COMB 3 successfully reproduces that the formation TiN is favored over the formation of Ti 2 N This is important for the atomic scale study of TiN alloy s Furthermore, the COMB3 Ti/TiN potentials correctly determine the work of adhesion on different interfacial structures of fcc Ti(001)/TiN(001) and to characterize the adsorption of oxygen atoms a nd molecules on the TiN(001) surface through deposition. In particular, the Ti/TiN potential could be easily integrated with other existing COMB
106 potentials, which facilitate variety studies of disparate systems at the large scaled MD simulations.
107 Table 5 1 COMB predictions compared to published experimental results, DFT calculations and MEAM potentials for the lattice constants a 0 and c 0 (), f H) (eV/atom), cohesive energy (E c ), bulk modulus (B), shear modulus (G) and elastic constants of the rocksalt (B1) TiN and rutile Ti 2 N (The reference states for the heats of formation are hcp Ti and gaseous N 2 ) Exp. DFT COMB 2NN MEAM b 2NN MEAM c TiN (B1) a 0 4.242 4.238 4.24 3 4.242 f H B1 1.74 1.70 1.7 5 1.74 1.82 E c 6.63 6.69 B 3 20 a 318 320 288 G 1 83 C 11 625 a 6 34 659 561 C 12 165 a 16 5 150 152 C 44 163 a 1 5 7 183 191 f H B4 f H B1 0.12 0.23 f H B3 f H B1 0.34 0.41 0.39 f H B2 f H B1 0.89 0.87 1.25 Ti 2 N (rutile) a 0 4.943 4.842 4.783 4.777 c 0 3.036 2.924 3.051 3.048 f H Rutile 1.38 1.26 1.26 1.63 1.63 f H Fe2N f H Rutile 0.72 0.23 0.015 f H Fluorite f H Rutile 0.62 0.34 f H Cu2O f H Rutile 0.68 0.45 Denotes fitted properties for COMB potential a Ref.[ 195 ] b Ref.[ 196 ] c Ref.[ 156 ]
108 Table 5 2 Comparison of formation energies (eV) of point defects in NaCl typed TiN crystal calculated from first principles calculations and COMB potentials. The chemical potentials of Ti and N are hcp Ti metal and an isolated N 2 molecule, respectively. Defects DFT COMB N vacancy 3.38 a 2.41 b,c 3.11 Ti vacancy 4.19 a 3.28 c 4.72 N interstitial 5.46 b 4.93 c 13.82 Ti interstitial 8.28 c 24.27 Schottky 5.28 5.80 c 7.22 Cation Frenkel 11.56 c 28.98 Anion Frenkel 7.34 c 16.93 a Ref.[ 158 ] b Ref.[ 197 ] c Ref.[ 198 ] Table 5 3 Relaxed surface energy of TiN (100), (110) and (111) surface predicted by COMB potentials, in comparison with the predictions from published first principles calculation and MEAM potentials. Surface (J/m 2 ) COMB MEAM First principles calculations (100) 1.41 1.30 a 1.66 b c (110) 2.07 2.46 a 3.01 b c (111) 2.54 3.65 a 5.28 b c a Ref.[ 199 ] b Ref.[ 156 ] c Ref.[ 157 ]
109 Figure 5 1. The unit cell of fitting structures for TiN and Ti 2 N systems in this study A ) NaCl B ) wurtzite C ) Zns D ) CsCl, E ) rutile, F ) Fe 2 N, G ) Cu 2 O, and H ) fluorite. Grey and blue atoms denote Ti and N atoms Figure 5 2. Schematic of layers sequence of the TiN(001), (110) (111) surfaces. Grey and blue atoms denote Ti and N atoms
110 A B Figure 5 3. The crystal structures of rocksalt TiN as a function of temperature. A) The atomic configurations where g rey and blue atoms denote Ti and N atoms B ) The pair distribution functions for the Ti N for rocksalt TiN from MD simulations at 300, 1000, 2000, 3000, and 3500 K.
111 A B Figure 5 4 The interfacial structures of Ti(001)/TiN(001), where Ti atoms in Ti metal are on the atop site of N, Ti and bridge site of Ti, respectively A ) DFT calculations; t he color scheme is the same as in Figure 1. B ) MD using COMB3; the atoms are color coded by atomic charge as indicated by the color bar
112 Figure 5 5 Snapshots of oxygen atoms de posited on the TiN(001) surface. A ) thirty and B ) sixty oxygen atoms deposited The color scheme is the same as in Figure 1. Figure 5 6. Snapshots from the deposition of single oxygen molecule on the TiN(001) surface. The atoms are color coded by the atomic charges with the charge values indicated by the color bar
113 CHAPTER 6 F IRST PRINCIPLES STUDY OF EARLY STAGE OXIDATION OF FACETED AlN S URFACES Alumina nitride (AlN) has great potential for use a thermal barrier material in protective coatings due to its properties of high theoretical thermal conductivity (320W/mK) and a small thermal expansion coefficient (4.6 x 10 6/ o C). It has been experimentally 200 202 demonstrated that when AlN is exposed to air at high temperature it is oxidized, influences its chemica l structure and thermal conductivity. Previous computational studies have investigated the adsorption of oxygen on flat AlN nonpolar 203 and polar 204 205 surfaces, but the influence of surface features and defects have yet to be taken into account in such ca lculations. For example, the reduced coordination of surface steps can substantially influence surface reactivity as has been demonstrated for the adsorption and dissociation of oxygen molecules on metal surfaces, 206 210 which show enhanced reactivity at steps in terms of smaller bond breaking energy barriers with respect to terrace sites. In this work, we use first principles, density functional theory (DFT) calculations to investigate the role of surface steps on the a dsorption of oxygen molecules and their subsequent dissociation on the AlN surface with a monoatomic step. To quantify the surface reactivity, the results are further compared to the energetics associated with a flat AlN surface. 6 1 Computational Details The DFT calculations are carried out using the Vienna ab initio simulation package (VASP). 59 62 The exchange correlation energy and potential are treated within the g eneralized gradient approximation (GGA) using the functional of Perdew Burke
114 Ernzerhof (PBE) 43 The valences electron orbitals for elements are Al3 s 2 3 p 1 N2 s 2 2 p 3 and O2 s 2 2 p 4 respectively. The energy cutoff for the plane wave basis is 500 eV. The AlN surface was built from the PBE optimized bulk AlN wurtzite structure. The properties of bulk AlN predicted by the DFT calculations are summarized in Table 6 1 and they compare well to published experimental values. The AlN surface structure is an eight layer slab with 64 atoms and is modeled using a 22 unit cell while t he corresponding stepped surface is modeled by a 43 system. A 14 vacuum space is added to separate the slabs, which i s determined to be sufficient to prevent sur face slabs from interacting with one another through the periodic boundary conditions The side and top views of the flat and stepped AlN surface structure are illustrated in Figure 6 1. In the calculations, the Brillouin zone i s s ampled with a 441 Monkhorst Pack (ref) k point mesh. The convergence criteria a re set at 1x10 5 eV and 1x10 2 eV 1 for energies and forces, respectively. For characterizing the interaction between oxygen molecule and AlN surfaces, t he adsorptio n energy is defined by the Equation 5 1, ( 6 1) where N is the number of O atoms adsorbed on the AlN surface, E O/AlN is the total energy of the oxyge n adsorbed on the surface, E AlN is the energy of the fully relaxed AlN surface, and E O 2 is the energy of single O 2 molecule in the vacuum. According to this equation, negative values mean bound states while positive values mean the states are unstable with respect to desorption.
115 6 2 Adsorption of Atomic and Molecular Oxygen on the AlN surface The preferred adoption sites are explored for atomic oxygen on the flat AlN surface by performing geometric optimization calculations. Three equilibrium adsorption sites denoted by A, B, and C in Figure 6 1 are identified for the (22) surface and Table 2 summaries the adsorption energ ies at thes e sites. Site A lies above the Al N dimer with an adsorption energy of 0.98 eV/atom, and site B lies between adjacent A N dimers with an adsorption energy of 0.59 eV/atom. Site C bonds two paralleled Al N dimers with adsorption energy of 0.10 eV/atom, w hich is much lower than the others It should be noted that a tomic oxygen placed right above the Al and N atom s are unstable and eventually move to the A site. In the case of atomic oxygen on the stepped AlN surface, based on the result of the preferred si te on the flat surface, three equilibrium adsorption sites that all lie above the Al N dimer and are denoted by D, E, and F as shown in Figure 6 1 are explored. The adsorption energies for these three sites are listed in Table 2 S ite F exhibit s the largest adsorption energy ( 1.08 eV/atom) followed by site D ( 1.00 eV/atom) and then site E ( 0.96 eV/atom) The results indicate that the step edge may lead to enhanced adsorption behavior. In the case of the adsorption and dissociation of the oxygen m olecule, the energy change for these processes is revealed in Figure 6 2. In Figure 6 2(a), the initial s tate starts with a flat surface with an O 2 lying far away from the surface (>8 ). When the O 2 is adsorbed on the surface, the total energy is decrease d by 1.37 eV/atom. The final state is O 2 dissociation and both oxygen atoms are sitting on the preferred site A, in which the total energy is decreased by 1.987 eV/atom compared to the initial state. As
116 for the stepped surface, the start of initial state is similar to that on flat surface, where O 2 l ies horizontally some distance from the surface. Because it is more complicated to anticipate the preferred sites for O 2 dissociation on a step ped surface, the oxygen molecule is placed above the Al N dimer in the vicinity of the top step edge to optimiz e the calculation, as illustrated in Figure 6 2(b). Following optimization, the results indicate that the oxygen molecule spontaneously dissociates into two oxygen atoms (d O O is 1.49 ) and the total energy is decreased by 2.096 eV/atom compared to the initial state. The spontaneous dissociation of O 2 is enhanced by the step edge where the undercoordinated surface atoms form bonds to the newly dissociated oxygen atoms Therefore, in the next step, we would like to discuss this spontaneous dissociation with electron density plots. 6 3 Summary In this work first principles, density functional theory calculations are used to investigate the role of surface steps on the adsorption of oxygen molecules and their subs equent dissociation The specific system considered is AlN with a monoatomic step. The results indicated that adsorption in the vicinity of the t op step edge leads to the spontaneous dissociation of O 2 which does not occur on the bare terrace. This work will serve as a benchmark for the testing and validation of future Al/Al 2 O 3 /AlN COMB potentials.
117 Table 6 1. Comparison of the computed bulk properties to the experimental data AlN (wurtzite) Experiment This work a () 3.110 3.109 ( 0.03 %) c () 4.978 4.990 (+0.24 %) a/c 1.601 1.608 (+0.43 %) E c (eV/atom) 7.497 7.446 ( 0.68 %) Table 6 2 Predicted adsorption energies (E ads )(eV) of atomic oxygen on the flat and stepped AlN surfaces. Adsorption site E ads ( flat ) E ads ( stepped ) top Al unstable top N unstable A 0.98 B 0.59 C 0.10 D 1.00 E 0.96 F 1.08
118 Figure 6 1. Stable adsorption sites on the flat and stepped AlN surface s. A) F lat surface structure and three stable adsorption sites (denoted by A B and C ). B) S tepped surface and three stable adsorption sites (denoted by D E and F )
119 A B Figure 6 2 The energy chage of the adsorption and dissociation of O 2 on the flat and stepped AlN surface A) flat and B) stepped.
120 CHAPTER 7 CONCLUSIONS In this dissertation, we report on the development and parameterization of third generation charge optimized many body (COMB) potential for the Cu/ZnO, TiN, and Ti/TiO 2 systems. Because the most important features of COMB are that it includes many body interactions with directional bonds and dynamic charge transfer, it is able to correctly model complex physical, chemical, and mechanical processes in these systems on their own or in combination with other systems parameterized within COMB. This approach was used to examine several different problems related to processes at surface s and interfaces. First, the growth mode of Cu on ZnO was characterized by modeling the deposition and subsequent Cu growth with COMB. Through MD simulations, the evolution of Cu growth was explored, and the results indicated that deposited Cu atoms form clusters on ZnO In particul ar, the simulation results clearly elucidated the ways in which Cu growth transitions from layer by layer to three dimensional as converge increases, a topic that is currently the subject of controversy in the literature. In addition, the influence of such factors as incident cluster deposition energy, ZnO surface temperature, and the nature of the oxide support on the potential nanoscaled changes at the interfaces such as Cu Zn mixing during high energy deposition and the different charge states of Cu ato ms on different supporting surfaces, were examined. These results are expected to provide useful guidance for the interpretation of experimental data and the optimization of supported metal systems for the catalyst applications.
121 Second, the properties of Ti, TiO 2 and Cu/TiO 2 were examined. The fitted and predicted properties of the Ti hcp metal and the TiO 2 rutile phase reasonably reproduced the experimental and density functional values. Most importantly, without the need for re parameterization but combi ng the Ti/TiO 2 potential developed in this work with a previously parameterized potential for Cu/Cu 2 O, COMB predicted that a preferred adsorption site for a single Cu atom is in between two bridging O atoms followed by the site above the bridging O atom an d the five fold coordinated Ti atom on the TiO 2 (110) surface, which agrees with the results of density functional theory calculations and experimental observations. Further, COMB predicted enhanced bonding between Cu clusters and the oxidized surface and i ndicated that a metal oxide bond between Cu and O was responsible for the high adsorption energy in this system. These results are expected to provide useful guidance to fabricate dispersed Cu clusters on TiO 2 and thus improve their catalytic activity. Fur thermore, the adsorption mechanisms identified in this work could be applicable to other transition metal/oxide systems. Third, the COMB potentials developed for TiN systems were used to investigate the relationship between interfacial structure and adhe sive properties of Ti/TiN interfaces. The results indicated that the COMB potentials correctly described the energetics of adhesion energy and charge distribution for three different interfacial structures of fcc Ti(001)/TiN(001), which are reasonably cons istent with density functional theory calculations. In addition, the characterization of the adsorption of oxygen atoms and molecules on the TiN(001) surface showed that the COMB potentials
122 for Ti/TiN/TiO 2 and N/NO systems were well suited to explore the c hemistry associated with their oxidation. Fourth, first principles, density functional theory calculations were used to investigate the role of surface steps on the adsorption of oxygen molecules and their subsequent dissociation on the AlN surface with a monoatomic step. The results indicated that adsorption at the vicinity of top step edge leads to a spontaneous dissociation of O 2 This work will serve as a benchmark for the testing and validation of Al/Al 2 O 3 /AlN COMB potentials c urrently being refined by other graduate students. While the COMB potential has a complex formalism for the potential energy, its design allows parameters to be flexible and transferable in multicomponent systems, as the examples in this dissertation illu strate. They are particularly well suited for the study of complex processes at surfaces and provide detailed atomic mechanisms that complement experimental observations and measurements. Most recently, the author participated in a more challenging task to examine the electrocatalytic reduction of CO 2 on the Cu/ZnO catalyst, as shown in Fig. 1 4, using the COMB potential. In this work, the products observed during MD simulations resemble those experimentally measured, which further demonstrates the ability of COMB to probe the complex process of surface catalysis. As a result of the work carried out and described in this dissertation, powerful new computational methods are thus available. These methods can be considered to be potent tools when used in combi nation with experimental methods to elucidate the details of complex processes in materials.
123 APPENDIX POTENTI A L PARAMETERS FOR COMB POTENTIALS DEVELOPED IN THIS WORK The atomic and electrostatic parameters of COMB potentials for pure elements (O, Cu, N, Ti, and Zn) are given in Table A 1 while t he bond typed parameters for pure elements and binary systems (Cu/ZnO, TiO 2 and TiN) are given in Tables A 2, A 3 and A 4
124 Table A 1. Atomic and electrostatic parameters of pure elements (O, Cu, N, Ti, and Zn) for COMB 3 potentials O Cu N Ti Zn X (eVq 1 ) 6.599630 3.652316 6.209731 3.095768 3.389925 J (eVq 2 ) 5.955097 3.213926 9.292255 4.230280 8.162758 K (eVq 3 ) 0.760433 0.747810 1.254466 1.039759 3.381140 L (eVq 4 ) 0.009388 0.117292 0.286350 0.357428 0.921681 1 ) 1.371794 1.397263 1.438711 0.724352 0.673439 Z (q) 1.53 9170 0.467753 0.552136 3.022932 0.620000 1 (eVq 3 r 3 ) 1.966411 0.309036 1.048322 0.365635 0.339814 1 ( eVq 4 r 5 ) 2.521788 0.926620 3.911805 0.764850 0.670207 D U () 1.213951 0.133947 0.688784 0.500000 0.50040 D L () 0.007664 0.349021 0.199693 0.005000 0.529610 Q U (q) 6 2 5 4 4 Q L (q) 2 2 3 4 4
125 Table A 2 Bond typed parameters of pure element (O, Cu, N, Ti, and Zn) for COMB 3 potential s O O Cu Cu N N Ti Ti Zn Zn A (eV) 4956.338821 712.352722 7654.972656 516.5873248 2842.164795 B 1 (eV) 688.163507 102.826096 2102.295654 117.0421345 494.811310 B 2 (eV) 0.000000 0.000000 0.000000 0.000000 510.776123 B 3 (eV) 0.000000 0.000000 0.000000 0.000000 530.021545 1 ) 5.295119 2.712035 5.218037 2.136011 2.626286 1 ( 1 ) 3.258854 1.467089 3.738549 1.178831 2.352653 2 ( 1 ) 0.000000 0.000000 0.000000 0.000000 2.330411 3 ( 1 ) 0.000000 0.000000 0.000000 0.000000 2.288269 3.258854 1.467089 3.738549 0.545629 2.364559 n 1.000000 1.000000 1.000000 0.566048 0.626739 m 1.000000 1.000000 1.000000 1.000000 1.000000 b 6 14.177827 0.000000 0.000000 0.550400 0.000000 b 5 16.001644 0.000000 0.000000 0.065549 0.000000 b 4 1.257097 0.000000 0.812665 0.130433 0.000000 b 3 5.652039 0.000000 1.724015 0.093054 0.000000 b 2 0.204688 0.000000 0.698596 0.226674 0.000000 b 1 1.826597 0.000000 0.283832 0.146606 0.000000 b 0 0.856557 0.231055 1.691325 0.077183 0.002774 c 0 0.000000 0.000000 0.000000 0.012318 0.002783 c 1 0.000000 0.000000 0.000000 0.175079 0.001013 c 2 0.000000 0.000000 0.000000 0.066085 0.002213 c 3 0.000000 0.000000 0.000000 0.076109 0.001191 n B 10.00 10.00 10 .00 10 .00 10.00 R s () 2.40 3.20 2.00 3.90 2.60 S s () 2.80 3.50 2.30 4.10 3.00
126 Table A 3 Bond typed parameters of Cu/ZnO system for COMB potential developed in this work. O Cu Cu O O Zn Zn O Cu Zn Zn Cu A (eV) 1336.465698 1336.465698 1407.226318 1407.226318 1318.298950 1318.298950 B 1 (eV) 934.717041 934.717041 155.372025 155.372025 416.695190 416.695190 B 2 (eV) 0.000000 0.000000 147.597641 147.597641 0.000000 0.000000 B 3 (eV) 0.000000 0.000000 177.190201 177.190201 0.000000 0.000000 1 ) 2.353223 2.353223 2.843915 2.843915 2.774643 2.774643 1 ( 1 ) 2.092751 2.092751 1.971994 1.971994 1.900813 1.900813 2 ( 1 ) 0.000000 0.000000 2.035631 2.035631 0.000000 0.000000 3 ( 1 ) 0.000000 0.000000 2.060415 2.060415 0.000000 0.000000 2.253174 2.229919 2.060761 0.725662 1.900813 1.900813 n 1.000000 1.000000 1.000000 0.626739 1.000000 0.626739 m 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 b 6 0.000000 0.002251 0.041682 0.024056 0.000000 0.000000 b 5 0.000000 0.001783 0.408499 0.002583 0.000000 0.000000 b 4 0.010825 0.018164 0.097052 0.067959 0.000000 0.000000 b 3 0.009698 0.003013 0.285974 0.028187 0.000000 0.000000 b 2 0.045080 0.002789 0.836569 0.067872 0.000000 0.000000 b 1 0.021738 0.004218 0.926345 0.110034 0.000000 0.000000 b 0 0.155773 0.000704 0.636860 0.050320 0.247469 0.227573 c 0 0.009366 0.034866 0.006728 0.036859 0.000000 0.000000 c 1 0.149421 0.012322 0.134903 0.047156 0.000000 0.000000 c 2 0.009016 0.019651 0.353478 0.000017 0.000000 0.000000 c 3 0.009550 0.021196 0.123085 0.015725 0.000000 0.000000 n B 10 .00 10 .00 10 .00 10 .00 10 .00 10 .00 R s () 2.30 2.30 2.60 2.60 3.30 3.30 S s () 2.60 2.60 3.00 3.00 3.6 0 3.6 0
127 Table A 4 Bond typed parameters of TiN and TiO 2 system s for COMB 3 potential s developed in this work. O N N O O Ti Ti O Ti N N Ti A (eV) 7654.972656 7654.972656 1556.205566 1556.205566 1876.810769 1876.810769 B 1 (eV) 2102.295654 2102.295654 212.525665 212.525665 194.957894 194.957894 B 2 (eV) 0.000000 0.000000 183.473007 183.473007 172.903008 172.903008 B 3 (eV) 0.000000 0.000000 181.943161 181.943161 198.120622 198.120622 1 ) 5.218037 5.218037 2.868830 2.868830 3.084135 3.084135 1 ( 1 ) 3.738549 3.738549 2.129506 2.129506 1.866843 1.866843 2 ( 1 ) 0.000000 0.000000 2.142030 2.142030 2.104133 2.104133 3 ( 1 ) 0.000000 0.000000 2.030774 2.030774 2.569293 2.569293 0.545629 3.738549 2.557314 0.791273 1.093141 3.364171 n 1.000000 1.000000 1.000000 0.566048 0.566048 1.000000 m 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 b 6 0.550400 0.000000 0.150761 0.049978 0.000002 0.000000 b 5 0.065549 0.000000 0.157712 0.154609 0.000022 0.000023 b 4 0.130433 0.812665 0.010905 0.044523 0.011337 0.021388 b 3 0.093054 1.724015 0.097477 0.169396 0.058494 0.266099 b 2 0.226674 0.698596 0.068680 0.179627 0.186848 0.539999 b 1 0.146606 0.283832 0.190148 0.239884 0.219402 0.103172 b 0 0.077183 1.691325 0.130799 0.184842 0.280179 0.061309 c 0 0.012318 0.000000 0.020522 0.003457 0.030106 0.052880 c 1 0.175079 0.000000 0.028074 0.022650 0.064553 0.034314 c 2 0.066085 0.000000 0.064522 0.016666 0.001688 0.020290 c 3 0.076109 0.000000 0.026573 0.003002 0.067116 0.020788 n B 10 .00 10 .00 10.00 10.00 10.00 10.00 R s () 2.00 2.00 2.80 2.80 2.60 2.60 S s () 2.30 2.30 3.20 3.20 3.00 3.00
128 REFERENCE LIST 1 U. Diebold, J. M. Pan, and T. E. Madey, Surf. Sci. 331 333 845 (1995). 2 V. E. Henrich and P. A. Cox, The surface science of metal oxides (Cambridge University Press, Cambridge, 1996). 3 H. O. Pierson ed., Handbook of refractory carbides and nitrides: properties, characteristics, processing, and applications (Noyes Publications, Park Ridge, N.J., 1996). 4 C. T. Campbell, Surf. Sci. Rep. 27 1 (1997). 5 U. Diebold, Surf. Sci. Rep. 48 53 (2003). 6 N. P. Padture, m. Gell, and E. H. Jordan, Science 296 280 (2002). 7 R. A. Miller, J. Am. Chem. Soc. 67 517 (1984). 8 M. Le, M. Ren, Z. Zhang, P. T. Sprunger, R. L. Kurtz, and J. C. Flake, J. Elect rochem. Soc. 158 E45 (2011). 9 P. Lazcano, M. Batzill, U. Diebold, and P. Hberle, Phys. Rev. B 77 035435 (2008). 10 G. Zhou and J. C. Yang, J. Mater. Res. 20 1684 (2005). 11 T. Do, S. J. Splinter, C. Chen, and N. S. Mclntyre, Surf. Sci. 387 192 (1997) 12 N. Birks, G. H. Meier, and F. S. Pettit, Introduction to the high temperature oxidation of metals (Cambridge University Press, Cambridge, 2009). 13 H. G. Karge and J. Weitkamp eds., Adsorption and diffusion (Springer, Verlag Berlin Heidelberg, 2008). 14 X. M. Liu, G. Q. Lu, Z. F. Yan, and J. Beltramini, Ind. Eng. Chem. Res. 42 6518 (2003). 15 T. Shishido, M. Yamamoto, D. Li, Y. Tian, H. Morioka, M. Honda, T. Sano, and K. Takehira, Appl. Catal., A 303 62 (2006).
129 16 Y. Tanaka, T. Utaka, R. Kikuchi, K. Sasaki, and K. Eguchi, Appl. Catal., A 238 11 (2003). 17 T. Shishido, Y. Yamamoto, H. Morioka, K. Takaki, and K. Takehira, Appl. Catal., A 263 249 (2004). 18 T. C. Germann and K. Kadau, Int. J. Mod. Phys. C 19 1315 (2 008). 19 J. E. Jones, Proc. Phys. Soc. 43 461 (1931). 20 J. E. Jones, Proc. R. Soc. Lond. A 106 463 (1924). 21 M. S. Daw and M. I. Baskes, Phys. Rev. Lett. 50 1285 (1983). 22 M. S. Daw and M. I. Baskes, Phys. Rev. B 29 6443 (1984). 23 R. A. Buckingham, Proc. R. Soc. London, Ser. A 168 264 (1938). 24 G. V. Lewis and C. R. A. Catlow, J. Phys. C: Solid State Phys. 18 1149 (1985). 25 M. W. Finnis and J. E. Sinclair, Phil. Mag. A 50 45 (1984). 26 F. Cleri and V. Rosato, Phys. Rev. B 48 22 (1993). 27 J. G Yu, S. B. Sinnott, and S. R. Phillpot, Phys. Rev. B 75 085311 (2007). 28 T. R. Shan, B. D. Devine, T. W. Kemper, S. B. Sinnott, and S. R. Phillpot, Physical Review B 81 125328 (2010). 29 B. Devine, T. R. Shan, Y. T. Cheng, A. J. H. McGaughey, M. Lee, S R. Phillpot, and S. B. Sinnott, Phys. Rev. B 84 17 (2011). 30 T. Liang, B. Devine, S. R. Phillpot, and S. B. Sinnott, J. Phys. Chem. A 116 7976 (2012). 31 A. C. T. van Duin, S. Dasgupta, F. Lorant, and W. A. Goddard, J. Phys. Chem. A 105 9396 (2001). 32 J. Tersoff, Phys. Rev. B 37 6991 (1988). 33 J. Tersoff, Phys. Rev. B 38 9902 (1988).
130 34 A. K. Rappe and W. A. Goddard III, J. Phys. Chem. 95 (1991). 35 T. R. Shan, B. D. Devine, S. R. Phillpot, and S. B. Sinnott, Physical Review B 83 115327 (2011). 3 6 D. W. Brenner, O. A. Shenderova, J. A. Harrison, S. J. Stuart, B. Ni, and S. B. Sinnott, J. Phys.: Condens. Matter 14 783 (2002). 37 K. H. Ernst, A. Ludviksson, R. Zhang, J. Yoshihara, and C. T. Campbell, Phys. Rev. B 47 13782 (1993). 38 J. M. Pomeroy, A. Couture, J. Jacobsen, C. C. Hill, J. P. Sethna, B. H. Cooper, and J. D. Brock, MRS Proc. 648 P7.3 (2000). 39 P. Hohenberg and W. Kohn, Phys. Rev. B 136 B864 (1964). 40 W. Kohn and L. J. Sham, Phys. Rev. 140 1133 (1965). 41 J. P. Perdew and A.Zunger, Phys. Rev. B 23 5048 (1981). 42 J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D.J. Singh, and C. Fiolhais, Phys. Rev. B 46 6671 (1992). 43 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77 3865 (1996). 44 D. Sholl and J. A. Steckel, Density Functional Theory: A Practical Introduction (Wiley New Jersey, 2009). 45 D. Frenkel and B. Smit, Understanding Molecular Simulation: From Algorithms to Applications (Academic Press, 2001). 46 T. Liang, B. Devine, S. R. Phillpot, and S. B. Sinnott, J. Phys. Chem. A 116 7976 (2012). 47 H. Toufar, K. Nulens, G. O. A. Janssens, W. J. Mortier, R. A. Schoonheydt, F. D. Proft, and P. Geerlings, J. Phys. Chem. 100 15383 (1996). 48 F. H. Streitz and J W. Mintmire, Phys. Rev. B 50 11996 (1994). 49 J. C. Slater, Phys. Rev. 36 57 (1930).
131 50 D. Wolf, P. Keblinski, S. R. Phillpot, and J. Eggebrecht, J. Chem. Phys. 110 8254 (1999). 51 P. P. Eward, Ann. Phys 369 253 (1921). 52 M. Wilson, M. Exner, Y. M. Huang, and M. W. Finnis, Phys. Rev. B 54 15683 (1996). 53 H. A. Stern, G. A. Kaminski, J. L. Banks, R. H. Zhou, B. J. Berne, and R. A. Friesner, J. Phys. Chem. B 103 4730 (1999). 54 J. Tersoff, Phys. Rev. Lett. 56 632 (1986). 55 A. Yasukawa, Jsme Int. J Ser. A 39 313 (1996). 56 D. W. Brenner, Phys. Rev. B 42 9458 (1990). 57 B. J. Thijsse, Nucl. Instr. and Meth. in Phys. Res. B 228 (2005). 58 J. A. Martinez, D. E. Yilmaz, T. Liang, S. B. Sinnott, and S. R. Phillpot, Curr. Opin. Solid State Mater. Sci., http://dx.doi.org/10.1016/j.cossms.2013.09.001 (2013). 59 G. Kresse and J. Hafner, Phys. Rev. B 47 558 (1993). 60 G. Kresse and J. Hafner, Phys. Rev. B 49 14251 (1994). 61 G. Kresse and J. Furthmller, Phys. Rev. B 54 11169 (1996). 62 G. Kresse and J. Furthmller, Comp. Mater. Sci. 6 15 (1996). 63 M. J. Frisch, et al., (Gaussian, Inc., Wallingford, CT, 2009). 64 W. E, W. Ren, and E. Vanden Eijnden, Phys. Rev. B 66 052301 (2002). 65 H. Jonsson, G. Mills, and K. W. Jacobsen, in Classical and Quantum Dynamics in Condensed Phase Simulations, Ed. B. J. Berne, G. Ciccotti and D. F. Coker, 385 (World Scientific, 1998). 66 G. Henkelman and H. Jonsson, J. Chem. Phys. 111 7010 (1999).
132 67 H. J. Ernst, F. Fabre, and J. Lapujoulade, Phys. Rev. B 46 1929 (1992) 68 M. Breeman and D. O. Boerma, Surf. Sci. Rep. 269 270 224 (1992). 69 L. Hansen, P. Stoltze, K. W. Jacobsen, and J. K. Nrskov, Phys. Rev. B 44 6523 (1991). 70 C. L. Liu, J. M. Cohen, J. B. Adams, and A. F. Voter, Surf. Sci. 253 334 (1991). 71 C. Lee G. T. Barkema, M. Breeman, A. Pasquarello, and R. Car, Surf. Sci. 306 L575 (1994). 72 G. A. Evangelakis and N. I. Papanicolaou, Surf. Sci. 347 376 (1996). 73 G. Henkelman and H. Jnsson, J. Chem. Phys. 111 7010 (1999). 74 S. Plimpton, J. Comp. Phys. 1 17 1 (1995). 75 A. O. Oluwajobi and X. Chen, Int. J. Auto. Comp. 8 326 (2011). 76 S. J. Plimpton and A. P. Thompson, MRS Bulletin 37 513 (2012). 77 G. C. Chinchen, K. C. Waugh, and D. A. Whan, Appl. Catal. 25 101 (1986). 78 J. Yoshihara and C. T. Campb ell, J. Catal. 161 776 (1996). 79 T. S. Askgaard, J. K. Nrskov, C. V. Ovesen, and P. Stoltze, J. Catal. 156 229 (1995). 80 M.M. Gnter, T. Ressler, B. Bems, C. Bscher, T. Genger, O. Hinrichsen, M. Muhler, and R. Schlgl, Catal. Lett. 71 37 (2001). 81 S. Sakahara, K. Yajima, R. Belosludov, S. Takami, M. Kubo, and A. Miyamoto, Appl. Surf. Sci. 189 253 (2002). 82 J. Nakamura, Y. Choi, and T. Fujitani, Top. Catal. 22 277 (2003). 83 J. B. Wagner, P. L. Hansen, A. M. Molenbroek, H. Topse, B.S.Clausen, and S. Helveg, J. Phys. Chem. B 107 7753 (2003).
133 84 J. D. Grunwaldt, A. M. Molenbroek, N. Y. Topse, T. Topse, and B. S. Clausen, J. Catal. 194 452 (2000). 85 T. Fujitani and J. Nakamura, Appl. Catal. A 191 111 (2000). 86 T. H. Fleisch and R. L. Mieville, J. Catal. 90 165 (1984). 87 J. K. W. Frese, J. Electrochem. Soc. 138 3338 (1991). 88 J. Flake (private communication). 89 O. Dulub, M. Batzill, and U. Diebold, Top. Catal. 36 65 (2005). 90 L. V. Koplitz, O. Dulub, and U Diebold, J. Phys. Chem. B 107 10583 (2003). 91 S. V. Didziulis, K. D. Butcher, S. L. Cohen, and E. I. Solomon, J. Am. Chem. Soc. 111 7110 (1989). 92 J. Yoshihara, J. M. Campbell, and C. T. Campbell, Surf. Sci. 406 235 (1998). 93 J. Yoshihara, S. C. Pa rker, and C. T. Campbell, Surf. Sci. 439 153 (1999). 94 O. Dulub, L. A. Boatner, and U. Diebold, Surf. Sci. 504 271 (2002). 95 G. Kresse, O. Dulub, and U. Diebold, Phys. Rev. B 68 245409 (2003). 96 A. Beltrn, J. Andrs, M. Calatayud, and J. B. L. Marti ns, Chem. Phys. Lett. 338 224 (2001). 97 D. R. Hummer, J. D. Kubicki, P. R. C. Kent, J. E. Post, and P. J. Heaney, Journal of Physical Chemistry C 113 4240 (2009). 98 K. Klier, Adv. Catal. 31 243 (1982). 99 A. C. T. v. Duin, S. Dasgupta, F. Lorant, and W. A. Goddard, J. Phys. Chem. A 105 9396 (2001). 100 D. Raymand, A. C. T. V. Duin, D. Spngberg, W. A. G. III, and K. Hermansson, Surf. Sci. 604 741 (2010).
134 101 D. Raymand, A. C. T. v. Duin, W. A. G. III, K. Hermansson, and D. Spngberg, J. Phys. Chem. C 115 8573 (2011). 102 A. C. T. v. Duin, V. S. Bryantsev, M. S. Diallo, W. A. Goddard, O. Rahaman, D. J. Doren, D. Raymand, and K. Hermansson, J. Phys. Chem. A 114 9507 (2010). 103 T. Liang, et al., Annu. Rev. Mater. Res. 43 109 (2013). 104 Y. K. Shin, T. R. Shan, T. Liang, M. J. Noordhoek, S. B. Sinnott, A. C. T. v. Duin, and S. R. Phillpot, MRS Bulletin 37 504 (2012). 105 V. F. Degtyareva, O. Degtyareva, M. K. Sakharov, N. I. Novokhatskaya, P. Dera, H. K. M ao, and R. J. Hemley, J. Phys.: Condens. Matter 17 7955 (2005). 106 Y. T. Cheng, T. R. Shan, B. Devine, D. Lee, T. Liang, B. B. Hinojosa, S. R. Phillpot, A. Asthagiri, and S. B. Sinnott, Surf. Sci. 606 (2012). 107 M. Hellstrm, D. Spngberg, K. Hermansson and P. Broqvist, Phys. Rev. B 86 235302 (2012). 108 G. Henkelman and H. Jnsson, J. Chem. Phys. 111 7010 (1999). 109 M. P. Allen and D. J. Tildesley eds., Computer simulation of liquids (Oxford univrsity press, New York, 1991). 110 L. Verlet, Phys. Rev 159 98 (1967). 111 J. Yoshihara, S. C. Parker, and C. T. Campbell, Surf. Sci. 439 153 (1999). 112 G. Ehrlich and F. G. Hudda, J. Chem. Phys. 44 1039 (1966). 113 R. L. Schwoebel and E. J. Shipsev, J. Appl. Phys. 37 3682 (1966). 114 C. R. Henry, Surf. Sci. Rep. 31 235 (1998). 115 K. H. Ernst, A. Ludviksson, R. Zhang, J. Yoshihara, and C. T. Campbell, Phys. Rev. B 47 13782 (1993). 116 O. Dulub, L. A. Boatner, and U. Diebold, Surf. Sci. 504 271 (2002). 117 L. V. Koplitz, O. Dulub, an d U. Diebold, J. Phys. Chem. B 107 10583 (2003).
135 118 U. Diebold, Surface Science Reports 48 53 (2003). 119 V. Swamy, J. D. Gale, and L. S. Dubrovinsky, Journal of Physics and Chemistry of Solids 62 887 (2001). 120 D. R. Collins, W. Smith, N. M. Harrison and T. R. Forester, Journal of Materials Chemistry 6 1385 (1996). 121 G. V. Lewis and C. R. A. Catlow, Journal of Physics C Solid State Physics 18 1149 (1985). 122 M. Matsui and M. Akaogi, Mol. Simul. 6 239 (1991). 123 R. T. Sanderson, Journal of the American Chemical Society 105 2259 (1983). 124 R. A. Donnelly and R. G. Parr, Journal of Chemical Physics 69 4431 (1978). 125 W. J. Mortier, K. Vangenechten, and J. Gasteiger, Journal of the American Chemical Society 107 829 (1985). 126 A. K. Rappe and W. A. Goddard, Journal of Physical Chemistry 95 3358 (1991). 127 A. Hallil, R. Tetot, F. Berthier, I. Braems, and J. Creuze, Physical Review B 73 165406 (2006). 128 R. Tetot, A. Hallil, J. Creuze, and I. Braems, Epl 83 (2008). 129 S. Y. Kim, N. Kumar, P. Persson, J. Sofo, A. C. T. van Duin, and J. D. Kubicki, Langmuir 29 7838 (2013). 130 S. Y. Kim and A. C. T. van Duin, The Journal of Physical Chemistry A 117 5655 (2013). 131 J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson M. R. Pederson, D. J. Singh, and C. Fiolhais, Physical Review B 46 6671 (1992). 132 H. J. Monkhorst and J. D. Pack, Physical Review B 13 5188 (1976). 133 G. Kresse and D. Joubert, Physical Review B 59 1758 (1999). 134 P. E. Blochl, Physical Review B 5 0 17953 (1994).
136 135 T. Liang, S. R. Phillpot, and S. B. Sinnott, Physical Review B 79 245110 (2009). 136 J. X. Zheng Johansson, O. Eriksson, and B. Johansson, Physical Review B 59 6131 (1999). 137 D. J. Bacon and V. Vitek, Metall. Mater. Trans. A 33 72 1 (2002). 138 J. D. Gale, Journal of the Chemical Society Faraday Transactions 93 629 (1997). 139 J. D. Gale and A. L. Rohl, Molecular Simulation 29 291 (2003). 140 E. G. Brovman, Y. Kagan, and A. Kholas, Soviet Physics Jetp Ussr 30 883 (1970). 141 M. Matsui, Journal of Chemical Physics 108 3304 (1998). 142 J. Scaranto, G. Mallia, S. Giorgianni, C. M. Zicovich Wilson, B. Civalleri, and N. M. Harrison, Surf. Sci. 600 305 (2006). 143 F. Labat, P. Baranek, and C. Adamo, J. Chem. Thoery Comput. 4 341 (2008). 144 C. R. A. Catlow, C. M. Freeman, and R. L. Royle, Physica 131B 1 (1985). 145 A. Hallil, E. Amzallag, S. Landron, and R. Tetot, Surface Science 605 738 (2011). 146 M. Hu, S. Noda, and H. Komiyama, Surf. Sci. 513 530 (2002). 147 B. Park, P. Liu, J. Hrbek, and J. F. Sanz, The Journal of Physical Chemistry C 113 7364 (2009). 148 L. Giordano, G. Pacchioni, T. Bredow, and J. F. Sanz, Surf. Sci. 471 21 (2001). 149 D. Phillay, Y. Wang, an d G. S. Hwang, Catal. Today 105 78 (2005). 150 Y. Tanizawa, T. Shido, W. J. Chun, K. Asakura, M. Nomura, and Y. Iwasawa, J. Phys. Chem. B 107 12917 (2003). 151 D. Matthey, J. G. Wang, S. Wendt, J. Matthiesen, R. Schaub, E. Laegsgaard, B. Hammer, and F. B esenbacher, Sciene 315 1692 (2007).
137 152 J. Hansen, et al., J. Phys. Chem. C 114 16964 (2010). 153 E. Lira, J. Hansen, L. R. Merte, P. T. Sprunger, Z. Li, F. Besenbacher, and S. Wendt, Top. Catal. 2013 1460 (2013). 154 R. P. Galhenage, H. Yan, S. A Tenney, N. Park, G. Henkelman, P. Albrecht, D. R. Mullins, and D. A. Chen, J. Phys. Chem. C 117 7191 (2013). 155 B. D. Devine, T. R. Shan, Y. T. Cheng, A. McGaughey, M. Y. Lee, S. R. Phillpot, and S. B. Sinnott, Phys. Rev. B 84 125308 (2011). 156 Y. M. Kim and B. J. Lee, Acta Mater. 56 3481 (2008). 157 H. Yu and F. Sun, Phys. B 404 1692 (2009). 158 Z. H. Xu, L. Yuan, D. B. Shan, and B. Guo, Comput. Mater. Sci. 50 1432 (2011). 159 D. G. Isaak, J. D. Carnes, O. L. Anderson, H. Cynn, and E. Hake, Physic s and Chemistry of Minerals 26 31 (1998). 160 H. Z. Yao, L. Z. Ouyang, and W. Y. Ching, Journal of the American Ceramic Society 90 3194 (2007). 161 J. Muscat, V. Swamy, and N. M. Harrison, Phys. Rev. B 65 224112 (2002). 162 P. E. Meagher and G. A. Lager, Can. Mineral 17 77 (1979). 163 W. Luo, S. F. Yang, Z. C. Wang, R. Ahuja, B. Johansson, J. Liu, and G. T. Zou, Solid State Commun. 133 49 (2005). 164 T. Mitsuhashi and O. J. Kleppa, J. Am. Ceram. Soc. 50 626 (1979). 165 L. E. Toth ed., Transition Metal Carbides and Nitrides (Academic Press, New York, 1971). 166 S. T. Oyama ed., The Chemistry of Transition Metal Carbides and Nitrides (Black Academic and Professional, London, UK, 1996). 167 P. Hedenqvist, M. Olsson, P. Wall n, Kassman, S. Hogmark, and S. Jacobson, Surf. Coat. Technol. 41 243 (1990). 168 H. E. Rebenne and D. G. Bhat, Surf. Coat. Technol. 63 1 (1994).
138 169 K. Hai, T. Sawase, H. Matsumura, M. Atsuta, K. Baba, and R. Hatada, J. Oral Rehabil. 27 361 (2000). 170 M. Wittmer, J. Vac. Sci. Technol. A 3 1797 (1985). 171 L. Hultman, Vac. 57 1 (2000). 172 S. B. Sinnott and D. W. Brenner, MRS Bulletin 37 469 (2012). 173 M. Marlo and V. Milman, Physical Review B 62 2899 (2000). 174 N. C. Hernndez and J. F. Sans, Comp. Mater. Sci. 35 183 (2006). 175 T. Iwasaki and H. Miura, J. Mater. Res. 16 1789 (2001). 176 Y. T. Cheng, T. R. Shan, T. Liang, R. K. Behera, S. R. Phillpot, and S. B. Sinnott, J. Phys.: Condens. Matter (submitted). 177 G. Kresse and J. Hafner, Physical Review B 47 558 (1993). 178 G. Kresse and J. Hafner, Physical Review B 49 14251 (1994). 179 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77 3865 (1996). 180 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 78 1396 (1997). 181 H. J. Monkhorst and J. D. Pack, Physical Review B 13 5188 (1976). 182 C. Lee, W. Yang, and R. G. Parr, Physical Review B 37 785 (1988). 183 A. D. Becke, J. Chem. Phys 98 5648 (1993). 184 M. Fuchs, J. L. F. D. Silva, C. St ampfl, J. Neugebauer, and M. Scheffler, Physical Review B 65 245212 (2002). 185 D. Wolf, Phys. Rev. Lett. 68 3315 (1992). 186 V. Fiorentini and M. Methfessel, J. Phys.: Condens. Matter 8 6525 (1996). 187 S. Plimpton, J. Comp. Phys. 117 1 (1995).
139 188 K. Harafuji, T. Tsuchiya, and K. Kawamura, J. Appl. Phys. 96 2501 (2004). 189 W. F. Wu, K. C. Tsai, C. G. Chao, J. C. Chen, and K. L. Ou, J. Electro. Mater. 34 1150 (2005). 190 L. A. S. Ries, D. S. Azambuja, and I. J. R. Baumvol, Surface and Coating Techno logy 89 114 (1997). 191 J. M. Lackner, L. Major, and M. Kot, Bulletin of the polish academy of science technical science 59 343 (2011). 192 J. Pelleg, L. Z. Zevin, and S. Lungo, Thin Solid Films 197 117 (1991). 193 S. V. Dudiy and B. I. Lundqvist, Physical Review B 69 125421 (2004). 194 J. Zimmermann, M. W. Finnis, and L. C. Ciacchi, J. Chem. Phys 130 134714 (2009). 195 R. G. Hennig, T. J. Lenosky, D. R. Trinkle, S. P. Rudin, and J. W. Wilkins, Physical Review B 78 054121 (2008). 196 J. O. Kim, J. D. Achenbach, P. B. Mirkarimi, M. Shinn, and S. A. Barnett, J. Appl. Phys. 72 1805 (1992). 197 S. S. Carara, L. A. Thesing, and P. Piquini, Thin Solid Films 515 2730 (2006). 198 L. Tsetseris, N. Kalfagiannis, S. Logothe tidis, and S. T. Pantelides, Physical Review B 76 224107 (2007). 199 R. Bs, Y. Pipon, N. Millard Pinard, S. Gavarini, and M. Freyss, Physical Review B 87 024104 (2013). 200 Y. Geng and M. G. Norton, J. Mater. Res. 14 2708 (1999). 201 J. W. Lee, I. Radu and M. Alexe, J. Mater. Sci. : Materials in electronics 13 131 (2002). 202 E. W. Osborne and M. G. Norton, J. Mater. Sci. 33 3859 (1998). 203 H. Ye, G. Chen, and Y. Wu, J. Phys. Chem. C 115 1882 (2011). 204 H. Ye, G. Chen, Y. Zhu, and S. H. Wei, Phys. Rev. B 77 033302 (2008).
140 205 M. S. Miao, P. G. Mosers, J. R. Weber, A. Janotti, and C. G. V. D. Walle, EPL 89 56004 (2010). 206 A. T. Gee and B. E. Hayden, J. Chem. Phys. 113 10333 (2000). 207 P. Gambardella, Z. Sljivancanin, B. Hammer, M. Blanc, K. Ku hnke, and K. Ken, Phys. Rev. Lett. 87 056103 (2001). 208 Y. Xu and M. Mavrikakis, Surf. Sci. 538 219 (2003). 209 Y. Xu and M. Mavrikakis, J. Phys. Chem. B 107 9298 (2003). 210 M. Ahonen, M. Hirsimki, A. Puisto, S. Auvinen, M. Valden, and M. Alatalo, Chem. Phys. Lett. 456 211 (2008).
141 BIOGRAPHICAL SKETCH Yu Ting Cheng was born in Taoyuan, Taiwan, in 1981 to Ming Shan Cheng and Li Hs iang Wu. He received his m echanical e ngineering and m aterials s cience and e ng ineering from National Sun Yat Sen University, Kaohsiung, Taiwan i n 2003 and 2005, respectively. In August 2008, Yu Ting was offered with the opportunity to pursue his Ph.D. degree at the University of Florida, and in May 2009 he was kindly provided with a research assistantship from Prof essor Susan B. Sinnott and since then working under her advisement towards the degree.