1 DEVELOPMENT AND APPLICATIONS OF VARIABLE CHARGE REACTIVE POTENTIALS FOR THE ATOMISTIC SIMULATIONS OF HETEROGENEOUS INTERFACES By TZU RAY SHAN 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 2011
2 2011 Tzu Ray Shan
3 To Sophia Peter and Maggie
4 ACKNOWLEDGMENTS Most dedicated appreciation to my advi sor, Prof. Susan Sinnott. I sincerely appreciate that I have had the opportunity to work with her. Her expertise in the field of computational materials science is truly admiring, and her passion for te aching and science has been my inspiration. Her endless patience and constant encouragement for me is really invaluable. I also deeply appreciate that I have had the honor to work with Prof. and Chair, Simon Phillpot. His class, Introduction to Molecular D ynamics, was definitely world class, and has opened up the world of simulations to me. During my PhD years, he worked almost side by side with Prof. Sinnott to instruct me and guide me through all the obstacles. They really are the best advisors one could ever wish for. This dissertation is also dedicated to my parents and my fianc e I would not have become who and what I am without the caring support from my parents. They backed me up in every decision I have made and supported me with anything I needed. The patience, blessing and encouragement of my fianc e Maggie, is what kept me going through these years. She is also my best friend and my soul mate. She took care of me in every way and enriched my life with happiness and laughter. What I have accomplis hed and achieved belongs to my parents and my fianc e Drs. Michael Chandross, Steve Plimpton and Aidan Thompson of the Sandia National Laboratories, New Mexico offered tremendous assistance during my three month visit in 2010. Dr. Chandross was a warm hos t during the visit and assisted in handling all kinds of tedious paper work. He also provided valuable comments and suggestions to my simulation works. Drs. Plimpton and Thompson assisted me in the coding of potentials; their knowledge in molecular dynamic s programs and empirical potentials is truly profound.
5 I would also like to thank my mentors from our g roup Dr. Wen Dung Hsu and Dr. Jianguo Yu. Wen Dung taught me the entire basics one should know of writing programs, running molecular dynamics simulations and parallel programming. I recall numerous after work hours that he has spent with me just to explain how things work. Jianguo, on the other hand lead me to the wonderful world (or endless abyss!) of potential development. He was the strictest mentor, yet the most brilliant and hardest working post doctoral researcher. Without them, my learning curve would be very slow and limited. Dr. Bryce Devin e has been the most important colleague throughout the course of study. Discussions and brain storms with him were what made all work possible. Along with Bryce, Drs. Tao Liang and Rakesh Behera are my go to guys whenever I have questions. I truly admire their deep knowledge in the fields of physics, chemistry, and computational materials science. I certainly enjoyed the time, both on and off work, spent with several members of the University of Florida Computational Materials Science Focus Group: Drs. Tak u Watanabe, Dilpuneet Aidhy, Donghwa Lee, Chan Woo Lee and Haixuan Xu, as well as Mr. Kun Ta Tsai, Travis Kemper and Eric Bucholz. Many thanks to the graduate, undergraduate, and high school students that I have mentored and worked with, Ya Ting Su, Karan Doshi, Scott Warren, Ira Hill, Yu Ting Cheng, Yanzhong Li, Jackelyn Martinez, Danielle Holton, Xuan Sun and Mark Noordhoek, for enriching my mentoring experience and even sharing my work load. Lastly, the funding source, National Science Foundation and Dep artment of Energy, and the computational resources managed by the Computational Materials Science Focus Group cluster administrators, are deeply appreciated.
6 TABLE OF CONTENTS page ACKNOWLEDGMENTS ................................ ................................ ................................ .. 4 LIST OF TABLES ................................ ................................ ................................ ............ 9 LIST OF FIGURES ................................ ................................ ................................ ........ 11 A BSTRACT ................................ ................................ ................................ ................... 13 CHAPTER 1 INTRODUCTION ................................ ................................ ................................ .... 15 2 MOLECULAR DYNAMICS, FUNCTIONAL FORMALISM AND IMPLEMENTATION OF CHARGE OPTIMIZED MANY BODY (COMB) POTENTIALS ................................ ................................ ................................ ......... 20 Fundamentals of Molecular Dynamics Simulations ................................ ................ 20 Inte r atomic Potential ................................ ................................ ....................... 21 Periodic Boundary Conditions ................................ ................................ .......... 21 Time Integration Algorithm ................................ ................................ ............... 22 Temper ature Control Method ................................ ................................ ........... 22 Pressure Control Method ................................ ................................ .................. 24 Functional Formalism of COMB Potentials ................................ ............................. 24 Self Energy Function ................................ ................................ ........................ 26 Inter atomic Potential Functions ................................ ................................ ....... 27 Charge Barrier Functions ................................ ................................ ................. 29 Correction Functions ................................ ................................ ........................ 30 Cosine bond bending function ................................ ................................ ... 30 Legendre polynomial bond bending function ................................ ............. 31 Additional short range repulsion ................................ ................................ 33 Over coordination function ................................ ................................ ......... 33 Tetragonal stabilizing function ................................ ................................ .... 34 Fitting Procedures of COMB Potentials ................................ ................................ .. 34 Dynamic Charge Equilibration Method ................................ ................................ .... 36 Implementation of COMB Potentials ................................ ................................ ....... 37 Implementation to In house Code ................................ ................................ .... 37 Implementation to Large scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) Code ................................ ................................ ........... 38 Scaling of LAMMPS Implementation of COMB Potentials ................................ 39 Scaling with size on single processor ................................ ........................ 39 Fixed size scaling on multiple processors ................................ .................. 40 Scaled size scaling on multiple processors ................................ ................ 41 Cost of COMB Potentials Compared to other Inter atomic Potentials in LAMMPS ................................ ................................ ................................ ....... 42
7 Summary ................................ ................................ ................................ ................ 43 3 DEVELOPMENT OF COMB SI/ SIO 2 POTENTIAL AND ITS APPLICATIONS TO AMORPHOUS SILICA AND SI/SIO 2 INTERFACES ................................ ............... 48 Potential Development and Properties of SiO 2 Polymorphs ................................ ... 49 Properties of Amorphous Silica and its Defect Formation Energies ........................ 51 Structural Properties of Amorphous Silica ................................ ........................ 52 Point Defect Formation Energies of Amorphous Silica ................................ ..... 52 Atomic Oxygen Adsorption Energies on Reconstructed Si (001) Surface .............. 53 Si/SiO 2 Interfaces ................................ ................................ ................................ .... 56 Quartz Interface ................................ ................................ ........................ 56 Si Nanocrystals Embedded in Amorphous Silica ................................ .............. 58 Nanoindentation of Si/SiO 2 Systems ................................ ................................ ....... 59 Summary ................................ ................................ ................................ ................ 62 4 DEVELOPMENT OF COMB CU/CU 2 O POTENTIAL AND ITS APPLICATION TO CU/SIO 2 INTERFACES ................................ ................................ .................... 75 Potential Development and Properties of Cu and Cu 2 O ................................ ......... 76 Adhesion of Cu/SiO 2 interfaces ................................ ................................ .............. 76 Ad cristobalite Interfaces ................................ ........................... 76 quartz Interfaces ................................ ................................ .. 80 Adhesion of Cu/a silica Interfaces ................................ ................................ .... 81 quartz Interfaces .............. 83 Summary ................................ ................................ ................................ ................ 84 5 DEVELOPMENT OF COMB HF/HFO 2 POTENTIAL AND ITS APPLICATIONS TO AMORPHOUS HfO 2 A ND SI/HFO 2 INTERFACES ................................ ........... 94 Potential Development and Properties of Hf Metal and HfO 2 Phases ..................... 95 Formalism and Fitting Procedures ................................ ................................ .... 95 Properties of Hf Metal ................................ ................................ ....................... 96 Properties of HfO 2 Phases ................................ ................................ ............... 98 Defect Formation Energies in Monoclinic HfO 2 ................................ .............. 101 Properties of Amorphous HfO 2 ................................ ................................ ....... 103 Applications to Si/HfO 2 Interfaces ................................ ................................ ......... 104 Si and Cubic HfO 2 Interface ................................ ................................ ........... 105 Structural properties of the interface ................................ ........................ 105 Adhesion of the interface ................................ ................................ ......... 106 Si and Monoclinic HfO 2 Interface ................................ ................................ .... 107 Nanoindentation of Si and Amorphous HfO 2 Systems ................................ .... 107 Summary ................................ ................................ ................................ .............. 108 6 DEVELOPMENT OF COMB TI /TIO 2 POTENTIAL AND ITS APPLICATIONS TO TIO 2 SURFACES AND TIO 2 SUPPORTED COPPER CLUSTERS AND ISLANDS ................................ ................................ ................................ .............. 120
8 Potential Development and Properties of Ti Metal and TiO 2 Phases .................... 121 Formalism and Fitting Procedure ................................ ................................ ... 121 Properties of Ti Metal ................................ ................................ ..................... 122 Properti es of TiO 2 Phases and Phase Stability ................................ .............. 123 Defect Formation Energies of Bulk Rutile TiO 2 ................................ ............... 125 Defect Formation Energies at Rutile TiO 2 Grain Boundaries .......................... 126 Surface Energies of TiO 2 Phases ................................ ................................ ......... 128 Low Miller Index Rutile Surfaces ................................ ................................ .... 128 Reconstructed Rutile (110) Surfaces ................................ .............................. 129 Step Rutile (110) Surfaces ................................ ................................ ............. 132 Low Miller Index Anatase Surfaces ................................ ................................ 134 Adsorption Energies of Cu n (n=2~4) Clusters on Rutile (110) Surface ................. 135 Adsorption Energies and Geometries ................................ ............................. 136 Potential Energy Surface of a Cu Dimer on Rutile (110) Surface ................... 137 Summary ................................ ................................ ................................ .............. 138 7 CONCLUSIONS ................................ ................................ ................................ ... 154 General Summaries ................................ ................................ .............................. 154 Ongoing Development of COMB Potentials ................................ .......................... 155 POTENTIAL PARAMETERS FOR COMB POTENTIALS DEVELOPED IN THIS W ORK ................................ ................................ ................................ ................... 157 REFERENCE LIST ................................ ................................ ................................ ...... 160 BIOGRAPHICAL SKETCH ................................ ................................ .......................... 173
9 LIST OF TABLES Table page 2 1 Relative computational cost of selected inter atomic potentials .......................... 44 3 1 cristobalite ............................... 63 3 2 quartz. ........... 64 3 3 Defect formation energies in amorphous silica ................................ ................... 64 3 4 Atomic oxygen adsorption energ ies on reconstructed Si ................................ .... 64 3 5 quartz and Si/a SiO 2 systems obtained from nanoindentation simulations ................................ ........................ 64 4 1 Properties of Cu 2 O ................................ ................................ ............................. 85 4 2 Adhesion energies of the fcc cristobalite (001) interfaces ................ 85 4 3 Adhesion energies of the fcc quartz (010) interfaces ....................... 85 4 4 Adhesion energies of the Cu/ a SiO 2 interfaces. ................................ ................. 85 4 5 Average Cu O and Cu Si bond lengths at the Cu/SiO 2 interfaces ...................... 86 5 1 Properties of hafnium metal ................................ ................................ .............. 109 5 2 Fitted properties of monoclinic, tetragonal and cubic hafnia phases ................ 110 5 3 Cohesive energy of monoclinic hafnia and that relative to monoclinic hafnia ... 111 5 4 Axial and volume thermal expansion coefficients of the monoclinic hafnia ....... 111 5 5 Defect formation energies of the monoclinic HfO 2 phase ................................ 111 5 6 Hardness and reduced modulus of Si/a HfO 2 structures. ................................ 111 6 1 Properties of Ti metal ................................ ................................ ....................... 139 6 2 Properties of rutile, anatase and brookite titania phases ................................ .. 140 6 3 Cohesive energies of rutile titania and that relative to rutile titania for other polymorphs ................................ ................................ ................................ ....... 142 6 4 Defect formation energies of bulk rutile. ................................ ........................... 142 6 5 Defect formation energies at rutile grain boundaries. ................................ ....... 143
10 6 6 Coordination number analysis of Ti and O atoms from the rutile structures containing Schottky or Frenkel defects ................................ ............................. 143 6 7 Surface energies of rutile surfaces. ................................ ................................ .. 143 6 8 Surface energies of various reconstructed rutile (110) surfaces. ...................... 143 6 9 Surface energies of rutile step (110) surfaces ................................ .................. 144 6 10 Surface energies of anatase surfaces. ................................ ............................. 144 6 11 Adsorption bond lengths and energies of Cu dimers, trimers, and tetramers ... 145 A 1 Potential parameters for COMB potentials developed in this work. .................. 158 A 2 Cutoff radii and mixing ru le scaling factors for short range interactions. ........... 159
11 LIST OF FIGURES Figure page 1 1 ................................ ................................ ............... 19 2 1 Scaling with size on single processor for LAMMPS Tersoff, LAMMPS COMB and HELL COMB ................................ ................................ ................................ 45 2 2 Fixed size scaling of COMB potential ................................ ................................ 46 2 3 Scaled size scaling of COMB potential ................................ ............................... 47 3 1 The cohesive energy as a function of unit volume curves for SiO 2 polymorph s 65 3 2 Amorphous silica obtained using COMB potential. ................................ ............. 66 3 3 The final relaxed structures of oxygen atoms adsorbed on a reconstructed Si (001) surface ................................ ................................ ................................ ...... 67 3 4 quartz interfacial structure after force and energy minimization ................................ ................................ ................................ ........ 68 3 5 quartz interfacial structure after equilibrating at 300K for 20ps ................................ ................................ ................................ ..... 69 3 6 quartz interface. ................................ ................................ ................................ ............. 70 3 7 Snapshot of a Si nanocrystal embedded in a S iO 2 matrix ................................ .. 71 3 8 Silanone bonds found at the Si/SiO 2 interface of the Si nanocrystal embedded in a SiO 2 matrix ................................ ................................ ................. 72 3 9 Snapshots of Si/a SiO 2 quartz systems from nanoindentation simulations ................................ ................................ ................................ .......... 73 3 10 Load displacement curves from nanoindentation simulations ............................ 74 4 1 Relaxed structures of the fcc cristobalite (001) interfaces ................ 87 4 2 Minimized structures of the fcc cristobalite (001) interfaces ............. 88 4 3 Minimized structures of an alternative intermediate oxygen concentration fcc cristobalite (001) interface ................................ ................................ 89 4 4 Minimized structures of the fcc quartz (010) interfaces .................... 90 4 5 Minimized structures of the fcc Cu (001)/a SiO 2 + 0 V O interfaces ..................... 91
12 4 6 Charge difference across the fcc cristobalite (001) interfaces .......... 92 4 7 Charge difference across the fcc quartz (010) interfaces. ................ 93 5 1 Cohesive energies as a function of unit volume for various hafnia phases ...... 112 5 2 Temperature dependence of the lattice constants of monoclinic hafnia ........... 113 5 3 Radial distribution function of amorphous hafnia ................................ .............. 114 5 4 Snapshots of the cubic HfO 2 and Si interfacial structure ................................ .. 115 5 5 Charge distribution of the c HfO 2 /Si interface for Si O and Hf ......................... 116 5 6 Snapshots of the moclinic HfO 2 and Si (m HfO 2 /Si) interfacial structure .......... 117 5 7 Snapshots of Si/a HfO 2 systems from nanoindentation simulations ................. 118 5 8 Load displacement curves of Si/a HfO 2 systems from nanoindentation simulations ................................ ................................ ................................ ........ 119 6 1 boundaries. ................................ ................................ ................................ ....... 146 6 2 Pristine and reconstructed rutile TiO 2 (110) surfaces ................................ ....... 148 6 3 Relaxed geometries of rutile (110) s tep surfaces ................................ ............. 150 6 4 Sphere representation of the (1x1) rutile TiO 2 (110) surface ............................ 151 6 5 Adsorption geometries of Cu on rutile surface ................................ .................. 152 6 6 Potential energy surface of a Cu dimer moving on (1x1) rutile (110) surface ... 153
13 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 DEVELOPMENT AND APPLICATIONS OF VARIABLE CHARGE REACTIVE POTENTIALS FOR THE ATOMISTIC SIMULATIONS OF HETEROGENEOUS INTERFACES By Tzu Ray Shan May 2011 Chair: Susan B. Sinnott Major: Materials Science and Engineering Electronic devices, gas sensors, micro and nano electromechanical systems and catalytic devices are usually composed of multifun ctional nano scaled structures. Inherent to t hese complex nanostructures are heterogeneous interfaces such as, for example, metal/metal oxide, semiconductor/metal oxide, and metal oxide/gaseous molecule interfaces. These heterogeneous interfaces play a crucial role in the applications of nanostructur es With recent advancements in technology, these nanostructures have shrunk to dimensions less than tens of nanometers and contain less than hundreds of billions of atoms. I t has also recently become routine to model and simulate these multi million atom nanostructures with heterogeneous interfaces at the atomistic level using molecular dynamics (MD) simulations The capability to simulate atoms interacting with each other in these multifunctional nanostructures has made probing, understanding and engineer ing the atomic level properties of materials possible. The key ingredient in MD simulations is the empirical energy functions used to des cribe interatomic interactions in such discrete bonding environments. E xamples of
14 successful empirical energy function s that can be used separately to describe interatomic interactions in metals, ionic solids and semiconductors are not flexible enough to model heterogeneous interfaces in multifunctional nanostructures. A r ecently developed first generation charge optimized many body (COMB) potential is designed to tackle this difficulty. In this work, the second generation COMB potential s are develop ed for various heterogeneous interfacial systems and are appl ied to investigate the underlying physics of the interfa ces in these nanostructures. The developed potentials include Si and SiO 2 Cu and Cu 2 O, Hf and HfO 2 and Ti and TiO 2 potential functions. The potential functions allow all elemental and oxide phases to be model ed in the same simulation cell, t herefore it i s able to model complex multifunctional nanostructures composed of heterogeneous interfaces in an integrated manner. The empirical, variable charge second generation COMB potentials developed in this work are useful and powerful tools that will allow scien tists and engineers to model real device size multifunctional nanostructures at the atomistic level with MD simulations using the COMB potentials, thus providing insights into the atomic level properties of materials with technological significance.
15 CHAPTER 1 INTRODUCTION Electronic devices, such as metal oxide semiconductor field effect transistors (MOSFETs), gas sensors, micro and nano electromechanical systems and catalytic devices are usually composed of multifunctional nano scaled structures. Inherent to these comple x nanostructures are heterogeneous interfaces such as, for example, metal/metal oxide, semiconductor/metal oxide, and metal oxide/gaseous molecule interfaces. It is typical for there to be significant changes in bonding environments on either side of thes e interfaces. For example, the bonding environment changes from metallic to ionic at metal/metal oxide interfaces, and changes from covalent to ionic at semiconductor/metal oxide interfaces. These heterogeneous interfaces play a crucial role in the applica tions of nanostructures: as illustrated in Figure 1 1 1 the application of thermal barrier coating materials involve interfaces between metallic and ionic bonded atoms, while interconnects and catalysts contain interfaces between metallic and covalent bonded atom s. In addition, interfaces between metallic, covalent and ionic bonding are common in microelectronic and bio implant applications. With recent advancements in technology, these nanostructures have shrunk to dimensions of several to tens of nanometers and contain only hundreds of millions to hundreds of billions of atoms. Driven by the rapid growth in computer power and computational algorithms, it has also recently become routine to model and simulate these multi million atom nanostructures with heterogen eous interfaces at the atomistic level. The capability to simulate atoms interacting with each other in these multifunctional nanostructures has made probing, understanding and engineering the atomic level properties of materials possible. The role of simu lations is not only
16 necessarily to model experimental behavior, but also to identify important atomistic mechanisms that occur during dynamical processes, to assist in the explanation and analysis of experimental data, and to provide predictions that conse quent experiments can confirm or refute. Reliable and dependable models and simulations allow scientists and engineers to minimize the amount of physical experimentation used for a study, reduce cost and avoid unnecessary generation of experimental waste. Modeling and simulating these heterogeneous interfaces in multifunctional nanostructured materials and devices with computers using ab initio quantum mechanical 2 and density functional theory (DFT) methods 3, 4 have provided valuable insights into the structural and energetic properties of the interfaces. The se types of electronic structure methods are effective and provide the highest fidelity currently available, but are also typically limited to systems with a few hundred atoms and are computationally time intensive. As a result of these limitations such ap proaches are not easily applied to the study of real device size d heterogeneous interfaces. Atomistic methods, most prominently, molecular dynamics (MD) simulations with empirical methods complement the strengths and weaknesses of electronic structure meth ods. MD simulations do not promise the same level of materials fidelity of electronic structure methods; however, they can simulate device size systems and dynami cal processes in a natural way. The key ingredient in MD simulations is the method of describi ng the way atoms interact with each other (or the empirical energy functions used to describe interatomic interactions) in such discrete bonding environments. The embedded atom method (EAM) 5, 6 fixed charge Buckingham type potentials 7 and Tersoff potentials 8, 9 are
17 examples of successful empirical energy functions that can be used to describe interatomic interactions in metals, ionic solids and semiconductors, respectively. However, to model and simulate heterogeneous interfaces in multifunctional nanostructures, the energy functions must be more flexible than they are in these potentials so that they are capable of describing all three types of bonding in an integrated manner. Furthermo re, the energy functions must allow for changes in atomic charge as a function of bonding environment, a seemingly inherently electronic effect. Reactive force fields (ReaxFF) 10, 11 and charge optimized man y body (COMB) 12, 13 potentials, therefore, are designed to ta ckle this difficulty and have been developed and parameterized to include a wide range of elements and compounds. In this regard, it is the objective of this work to develop and apply COMB potentials for a various heterogeneous interfacial systems and inve stigate the underlying physics of the interfaces in these nanostructures. The rest of the chapters of the dissertation are organized as follows. Chapter 2 reviews the fundamentals of MD simulations, the details of the inter atomic empirical potential funct ion used in this dissertation, dynamic charge equilibration methodology and the implementation of the potential into MD codes. Chapter 3 discusses a potential function for semiconductor Si and gate oxide SiO 2 polymorphs. The potential function for Si/SiO 2 is applied to investigate the point defect formation energies of amorphous silica, atomic oxygen adsorption energies on reconstructed Si (001) surface, Si and SiO 2 interfaces, Si nanocrystals embedded in amorphous silica matrix, and nanoindentation of Si/ quartz and Si/a SiO 2 systems. Chapter 4 provides a potential function for metal Cu and Cu 2 O oxide. The potential for Cu/Cu 2 O, when combined with
18 the potential for Si/SiO 2 is applied to study the structural and adhesive properties of Cu/SiO 2 interfaces. Chapter 5 gives a potential function for metal Hf and gate oxide HfO 2 polymorphs. The Hf/HfO 2 potential captures the correct phase order among HfO 2 polymorphs and is applied to study amorphous hafnia. The potential for HfO 2 together with Si/SiO 2 potential is also applied to investigate the structural and adhesive properties of Si/HfO 2 interfaces, as well as nanoindentation of Si/HfO 2 interfacial systems. Chapter 6 is a potential function for metal Ti and metal oxide TiO 2 polymorphs, which is applied to in vestigate the phase order of TiO 2 polymorphs, defect formation energies in bulk rutile and in rutile grain boundaries, surface energies of various rutile and anatase surfaces. The Ti/TiO 2 potential is also combined with the Cu/Cu 2 O potential to study the a dsorption energetic of Cu clusters on rutile (110) surfaces. Lastly, Chapter 7 contains a general conclusion of this work, along with ongoing development of new potential functions and future possibilities.
19 Figure 1 1. described in a recently published article ( Science 325, 1634 (2009) ) illustrates the significant difference in bonding environments between interfaces of multifunctional nanostructures Reprinted with permission from S. R. Phillpot and S. B. Sinnott
20 CHAPTER 2 MOLECULAR DYNAMICS, FUNCTIONAL FORMALISM AND IMPLEMENTATION OF CHARGE OPTIMIZED MAN Y BODY (COMB) POTENT IALS 2.1 Fundamentals of Molecular Dynamics Simulations Molecular dynamics is a computer simulation technique that the evolution of a set of in teracting atoms with respect to time is followed by integrating their equation s of motion. 14 equation of motion is numerically integrated to calculate atomic posit ions and momenta: (2 1) (2 2) where is the atomic mass of atom i is its acceleration, is its position and is the force acting upon it, due to inter atomic interactions. The bold text represents a vector, and this representation will be used throughout this chapter. MD is a statistical mechanics method, and it provides a way to obtain a se t of configurations (atomic positions and momenta) distributed according to some statistical distribution function (ensemble). From statistical mechanics 14 the average over such a configuration distributed according to a certain ense mble represents a physical quantity. Statistical mechanics also serves as a link between microscopic behavior and macroscopic thermodynamics. Therefore provided a long enough simulation time for a system of atoms to reach equilibrium, MD simulations can be used to measure thermodynamic properties of the system and the material simulated. For further details on MD simulations and statistical mechanics, the reader is referred to Allen [ 14 ]. Only the basics of MD simulations, s uch as inter atomic potential, periodic boundary conditions,
21 time integration algorithm and temperature/pressure control methods, are briefly introduced in the following sections. 2.1.1 Inter atomic Potential MD simulations utilize energy functions to describe inter atomic interactions, and such energy functions, or the inter atomic potential, governs the force on atom i by: (2 3) where is the inter atomic potential energy function. The inter atomic potential can be as simple as a pairwise interaction between a pair of atoms, which is the case in a Lennard Jones potential 15 or includes bond angles between three bonded atoms, or even dihedral angles between four bonded atoms. The potential function may also include long range interactions between non bonded atoms. The potential may depend only on atomic positions or it may be also dependent upon atomic charges, q i The choice of inter atomic potential is crucial, since the quality of MD simulations and the fidelity of the thermodynamics properties obtained for the system/material rely heavily on the in ter atomic potential. The inter atomic potential used in this study is the charge optimized many body 12, 16, 17 (COMB) potential, the det ails of which will be discussed in Section 2.2. 2.1.2 Periodic Boundary Conditions Periodic boundary conditions (PBCs) 14 are introduced to avoid unphysical edge effects, and allow a finite number of atoms, distributed in a simulation supercell, to model an infinite system. In this case, the simulation supercell is surrounded by replicas of itself so that atoms that exit one side of the supercell remerge into the simulation through the opposite side of the supercell. Moreover, an atom in the supercell interacts
22 with either another atom in the supercell or its equivalent atom in a surrounding imaginary supercell depending on which distance to the atom is shortest. I t should be noted that PBCs are approximations, and the influence of boundaries on the system is not completely eliminated. However, the undesired influence of boundaries can be minimized by choosing sufficiently large dimensions for the simulated system. 2.1.3 Time Integration Algorithm The evolution of a system containing N atoms with time t in MD simulations occurs by integrating Equation 2 1 with a time integration algorithm. In particular, at a certain time t knowing the positions of the atom i the positions of the atom at a later time t t can be obtained. Through iteration the evolution of the system of N atoms can be followed for long times. In this work the time integration algorithm utilized is the Verlet algorithm. 14, 18 The algorithm expands the positions, to two third order Taylor expansions in to obtain the positions at a later time: (2 4) where a equation of motion is simply: (2 5) The tru ncation error of this algorithm is of the order of With a carefully chosen time step, that is small enough the truncation error can be neglected. 2.1.4 Temperature Control Method Thermostats are often used i n MD simulations to regulate system temperature. Since temperature is an averaged physical quantity, one cannot control the system temperature directly by fixing it to a certain value. The system temperature is therefore controlled by regulating the kineti c energy of the system; in practice, one modifies the
23 velocities of the atoms. Of the several popular thermostat methods 19 23 the Nos Hoover 20 22 thermostat is applied in this work. In the Nos Hoover thermostat formalism, the effect that modifies the atomic velocities, which acts as a heat reservoir, is reduced equation of motion, by the following equation s: (2 6) (2 7) Where p is the momenta of the atoms, and is the additional degree of freedom representing the heat reservoir. The equation of motion for the reservoir is: (2 8) where is the number of independent momentum degrees of freedom of the system, is the Boltzmann constant, T is the system temperature and T 0 is the reservoir temperature. Q strength of the thermostat, and should be chosen carefully. If the value of Q is too small, the system will not behave as a canonical, i.e. constant volume and constant temperature (NVT), e nsemble. However if Q is too large then the temperature control will not be efficient. In MD simulations, one can selectively apply thermostat to atoms of the system. Atoms that are not subject to the constraint of the thermostat, termed active atoms, are allowed to evolve solely dependent upon inter atomic potential. Selectively applying thermostat allows dynamic properties and chemical reactions being captured, but system temperature is efficiently regulated by the reservoir. The Nos Hoover thermostat is applied in this work for MD simulations carried out in NVT ensemble.
24 2.1.5 Pressure Control Method Two of the most considered ensembles, canonical and microcanonical, in MD simulations are carried out under constant volume constraints. In such simulations the volume of the system containing atoms remains constant. However for some cases one would wish the volume of the system respond to external effects such as temperature, in which case one controls the pressure of the system and allows the volume to cha nge. These pressure control methods are called barostats, and several barostats have been developed for MD simulations, as seen in the literature. 19, 23, 24 This work utilizes Martyna Tobias Klein (MTK) implementation 24 of Nos Hoover barostat, which has the form: (2 9) (2 10) where is the additional d change of volume, P is the system pressure, P 0 is the piston pressure and p is the barostat time constant that controls the strength of the pressure control. Because of its implementation, MTK Nos Hoover barostat has to be used together with the Nos Hoover thermostat, which forms the constant pressure constant volume (NPT) ensemble. This combination of Nos Hoover thermostat and MTK Nos Hoover barostat is applied in this work for MD simulatio ns carried out in NPT ensemble. 2. 2 Functional F ormalism of COMB P otentials The COMB potential s developed in this work are built on the physical and mathematical framework developed by others for the Si/SiO 2 system. 12 In brief the first generation COMB potential formalism is based on the empirical Tersoff potential for
25 silicon. 8, 9 The Tersoff potential describes many body interactions, allows the breaking of existing bonds and the formation of new bonds, and has been successfully applied to study various properties of silicon and diamond. B ased on the Tersoff potential, Yasukawa extended the bond order method to the Si/SiO 2 system 25 by adding an electrostatic term and self consistent charge equilibration in the spirit of the Rappe and Goddard approach 26 which was also applied to Al/Al 2 O 3 by Streitz and Mintmire 27 Iwasaki applied the same Yasukawa approach to the Si/HfO 2 Ge/HfO 2 and Si/ZrO 2 systems, 28, 29 while Yu et al. modified the Yasukawa potential function by adding correction terms, yielding the first generation COMB framework for Si/SiO 2 12 Although these extended Tersoff (Yasukawa, Iwasaki and Yu) potentials have their strengths, including the ability to model multi component heterogeneous systems and transferability among different species, they are limited because of the use of a damped point charg e model and because of the unrealistically short cutoff for the electrostatic term in the potential. To rectify these deficiencies, the author of this work and others modified the potential function for the Cu/Cu 2 O systems 30 by replacing the point charge model and the cutoff function with Coulomb integrals over Slater 1s orbital 31 and by treating the electrostatic interactions using the real space Wolf summation method. 32 Similar to Brenner deriving the R eactive Empirical Bond Order (REBO) 33, 34 potential from DFT 35 and Finnis deriving an effective n body potential from t ight binding DFT, 36 Devine show ed that the second generation COMB potentials can be derived from DFT as well. 17 This work, bases on the fundamental s described above, develops the second generation COMB potenti al formalism for Si/SiO 2 16 Cu/Cu 2 O 37 Hf/HfO 2 13 and Ti/TiO 2 38 The second generation COMB potential s ha ve the general functional form:
26 ( 2 11 ) where is the total potential energy of the system, is the self energy term of atom is the inter atomic potential between the th and th atoms, is the distance between atoms and and are charges of the atoms and and and are the charge barrier and correction functions The individual functions are further described in the following sections. 2.2. 1 Self Energy Function The self energy term describes the energy required to form a charge and takes the form: ( 2 1 2 ) where the coefficients and are adjustable parameters fit to the atomic ionization energies and electron affinities of oxygen, silicon, copper, hafnium and titanium The oxygen coefficients are the same for all COMB potentials (Si/SiO 2 16 Hf/HfO 2 13 Cu/Cu 2 O 37 and Ti/TiO 2 38 ) to ensure full compatibility and transferability of the different potentials. A penalty function that captures the change in atomic self energy due to the field of the ionic lattice such as bulk oxides, is added to the self energy term and takes the form: ( 2 1 3 ) where and are adjustable parameters fit to the bulk properties
27 2.2. 2 Inter atomic Potential Functions The inter atomic potential energy consists of three components: short range repulsion, short range attraction, and long range Coulombic interaction, which are defined as: (2 1 4 ) (2 1 5 ) (2 1 6 ) (2 1 7 ) The many body effects are described with the bond order term, in the short range attraction, which takes the form ( 2 1 8 ) where the symmetry function and angular function are defined as ( 2 1 9 ) (2 2 0 ) Here is the bond angle between bonds and , and are adjustable parameters. The inverse decay lengths and the energy coefficients and depend only on the types of interacting atoms and follow take the following forms:
28 ( 2 2 1 ) ( 2 2 2 ) ( 2 2 3 ) ( 2 2 4 ) where and are adjustable mixing rule scaling factors. A change in the partial charge on an atom affects the effective ionic radius, which influences the short range attraction and repulsion. Yasukawa 25 addresse d this issue by modifying and to be depend ent on the charge of the atom : ( 2 2 5 ) ( 2 2 6 ) ( 2 2 7 ) ( 2 2 8 ) ( 2 2 9 ) ( 2 3 0 ) ( 2 3 1 ) ( 2 3 2 ) ( 2 3 3 )
29 where and are adjustable parameters that reflect the change in ionic radius with charge, and are parameters that correspond to upper and lower limits of valence shells and is an adjustable parameter. The short range repulsion and attraction are truncated by the cut off function which is defined as: ( 2 3 4 ) where R and S are optimized cutoff radii that truncates after first nearest neighbors. The long range Coulombic interaction between charged atoms is described with the charg e coupling factor, and takes the form : ( 2 3 5 ) ( 2 3 6 ) which is a Coulomb integral over 1 s type Slater orbitals 31 ; is an orbital exponent that controls the radial decay of the spherical electron density. The Coulombic interaction is treated with the Wolf summation 32 and takes a cutoff of 12 .0 2.2. 3 Charge Barrier Functions The charge equilibration method although solves atomic charges dynamically, assumes that all atomic charges are limited to the valence shell, which are defined by the parameters and Once the atomic charge during charge equilibration reaches the upper or lower limit of the valence shell, the charge dependent energy
30 functions exhibit an abrupt discontinuity within/without the shells (due to Equation s 2 2 5 to 2 3 3) Since the potential function requires continuous c harge dependent terms c harge barrier functions are applied when charges reach the valence shell limits and take the form: (2 3 7) (2 3 8) w here and are adjustable energy parameters currently set to 1keV Different from the self energy function and inter atomic potential functions, the charge barrier function is not derived from DFT. The choice of a fourth order polynomial for the charge barrier func tion is arbitrary. However, for aesthetic reason the fourth order polynomial was chosen to match the fourth order self energy term (described in Section 2.2.1, Equation 2 1 2). 2.2. 4 Correction Functions Unlike the self energy, inter atomic potential, and charge barrier functions discussed in previous sections correction functions are not applied to each individual potential as a general practice. Ideally speaking, correction functions should be avoided since they are not directly related to chemical bondi ng and they may not have significant physical meaning. However, these functions are applied only when it was found to be necessary to yield a more robust potential These corrections are applied to bonded pairs and bond angles as described in the following sections 2.2. 4 .1 Cosine bond bending function The cosine function bond bending terms are generally applied to stabilize the ground state phase relative to other less er stable polymorphs. The function is only
31 applied to oxygen containing bonds with un like pairs. Although the angular function has a similar purpose, it ( Equation 2 3 0) is applied to all bonds including M M M, M O M, O M O, O O M, and O O O bonds, where M is the metal/ cation specie Therefore the cosine bond bending terms, applied solely to M O M and O M O bonds, are defined as: ( 2 3 9 ) ( 2 4 0 ) where M are Si, Cu and Ti, and are adjustable bond angle strengths, and and are set to the ideal bond angles of the ground state phase s of SiO 2 Cu 2 O and TiO 2 This function is not applied to the COMB Hf/ HfO 2 potential for the complexity of Hf O Hf and O Hf O bond angles, which will be further discussed in Section 2.2. 4.4. 2.2. 4 .2 Legendre p olynomial bond bending function The Legendre polynomial bond bending function applied to M M M bonds, are also used to stabilize ground state phases. Since the short range interactions are truncated after first nearest neighbors and the Coulombic contribution is negligible because of neutral charges for pure metal phases, the COMB potential function do es not distinguish between a metal in the h exagonal close packed ( hcp ) or face centered cubic ( fcc ) phases without strong coefficients for the angular term, Equation 2 3 0. However increased strengths of the angular term to capture the energy difference between hcp and fcc phases compromises the mec hanical properties. The hcp phase is only significantly different from the fcc phase when third nearest neighbors are
32 considered. However inclusion of explicit third nearest neighbor interactions within an empirical potential significantly increases its computational cost. To cope with this difficulty, Yu et al. applied a third order Legendre polynomial, that was first proposed by Th ijsse 39 to Cu Cu Cu bonds in pure Cu 40 which takes the form: ( 2 4 1) where is the third order Legendre Polynomial function of the Cu Cu Cu bond angle, and is the coefficient fitted to the difference in cohesive energies of Cu in hcp and fcc phases. This function penalizes the 109.47 angle of the hcp phase in order to make the fcc phase more stable and is used only in the COMB Cu/Cu 2 O potential so far It is expected for future COMB potentials of fcc metals to include this bond angle correction function. When applying this function to hcp metals such as Ti and Hf, however, the leading coefficient must be negative since hcp is the more stable phase. Also, the absolute value of the fcc hcp energy difference for Ti/ Hf is approximately an order of magnitude greater than that of Cu; thus must be a larger negative number. Although this function with a negative leading coefficient does successfully stabilize the 109.47 bond angle of the hcp phase, it also penalizes the 146 angle that is specific to the hcp phase. This penalty is significant because of the large value, which leads to an unstable hcp phase above approximately 200 K. To overcome this problem while still retai ning a more stable hcp phase, Shan et al. in their COMB Hf / Hf O 2 13 and Ti / Ti O 2 38 potentials employ ed a sixth order Legendre p olynomial to metal bonds which is defined as:
33 (2 4 2) where M are Ti and Hf, is the sixth order Legendre polynomial function of the bond angle, and is the coefficient fitted to the difference in cohesive energies of Ti/ Hf in hcp and fcc phases. This function stabilizes the 146 bond angle of the hcp phase, wi thout significantly penalizing the 109.47 angle. With this polynomial function, the hcp phase of Hf is stable and the hcp fcc energy difference is well described. This sixth order Legendre polynomial bond angle correction function is also expected to be a pplied to future COMB potentials of hcp metals. 2.2. 4 .3 Additional short range repulsion An additional short range repulsion term, applied only to Ti O bonds to bring the equilibrium bond lengths and cohesive energies of TiO 2 phases closer to first principles calculated values, is added to the short range repulsion and is defined as: ( 2 4 3) where is the fitted strength and is set to S the outer short range cutoff radius. 2.2. 4.4 Over coordination function The challenge for the HfO 2 system is the variety of Hf O Hf and O Hf O bond angles that are displayed by the various hafnia phases. For this reason, attempts to stabilize the ground phase with bond bending terms applied to Hf O Hf and O Hf O at first proved unsuccessful. Therefore, I introduce an over coordination correction term to de stabilize the phases with higher (8 and greater) coordination number around the Hf atom. This term was inspired by the analogous term in the ReaxFF potential for hydrocarbons 10 and takes the form for the Hf O system:
34 ( 2 4 4 ) ( 2 4 5 ) Here E 0 and are fitted parameters that control the strength of the over coordination correction, and is the over coordination number with respect to 7 (the coordination number of Hf atom in the ground state monoclinic phase). The calculation of incorporates the short range cutoff defined in Equation 2 3 4 where R S and S S are th e inner and outer covalent cutoffs, respectively. This correction term is only applied when is greater than 0.25. 2.2. 4 5 Tetragonal stabilizing function To stabilize the tetragonal HfO 2 phase over the lesser stable cubic phase, an additional repulsion term that takes the following form is applied to Hf O bonds : ( 2 4 6) ( 2 4 7 ) where t he coefficient s r 1 and r 2 are adjustable parameters and r Hf O is the Hf O bond length. This additional term ( Equation 2 4 6) is added into the repulsive energy term ( 2 1 5 ). 2. 3 Fitting Procedures of COMB Potentials The parameter sets for the previously developed first generation COMB potentials for Si/SiO 2 12 and Cu 40 were fit s olely to the structural and energetic information of the ground state polymorphs. Yu et al. determined the best parameter set that minimize a penalty function, f(p) for the parameter set, p with the least squares method using a downhill simplex algorithm 41 The penalty function is calculated from the weighted
35 deviations of the lattice parameters, elastic constants, and energy per isometric strain of the ground state phase from respective target values obtained from experiments or first principles calculations. These properties of all other less stable polymorphs are then checked against the resulting parameter set. Since the param eters are only optimized for the ground state phase, the drawback of this method is that the optimized parameter set must be re determined if some less stable polymorph is incorrectly predicted to be lower in energy than the ground state or exhibits unacc eptable properties. To overcome this difficulty, the author of this work ha s developed an algorithm that fits structural and energetic information for multiple phases. It also determines the potential parameter set using the least squares method with the s implex algorithm, but t he penalty function is det ermined from all considered phases which generally include natural polymorphs, high temperature/pressure polymorphs, and hypothetical stoichiometric and non stoichiometric polymorphs. It is also possible an d feasible to include structures with defects, both point defects and planar defects such as free surfaces and stacking faults. A major advantage of this multi phase fitting method compared to the ground state method is that the fitting of parameters consi ders all possible structures and configurations including perfect and defective ones so that the relative difference in cohesive energies and the phase order is at least taken into account, if not well reproduced. This method has been and still is applie d to the development of all second generation COMB potentials. 13, 16, 30, 37, 38, 42 45 The disadvantage, as one might expect is the increase in computational expense within each simplex fitting iteration The parameters determined for the potentials developed or modified as part of this work are given in Appendix, Tables A 1 and A 2.
36 2. 4 Dynamic Charge Equilibration Method Atomic charges are treated as dynamica l variables that evolve with time and bonding e nvi ronment in COMB potentials, and are equilibrated (solved) based on the electronegativity equalization (EE) principle a concept after Sanderson 46 and Parr 47, 48 Sanderson proposed that in a many atom system, the electron s, regarded as an electron gas will distribute themselves so that the electrochemical potential is equal at all atomic/nuclear sites. It wa s later shown by Parr et al. that the electron gas has a n electrochemical potential which is the negative of the Mulliken electronegativity. The atomic charges on all atoms, by the EE principle, are those that minimize the charge dependent energy functions and for which the electronegativities are equal (or within an arbitrary tolerance factor). If an atom/nucleus moves, the electron gas will re equilibrate itself to equalize the electronegativi ty on all atoms which results in a change in atomic charges. Treating charges as dynamical variables adds an additional degree of freedom to the atoms besides the positions. The main task for the potential is to explicitly solve the positions and charges correctly and efficiently. Traditionally, for an N atom system a self consistent iterative looping me thod involving an N N matrix is used to solve the charges. This method, however, also scales with O(N N), which can significantly increase computational cost. To bypass this problem, the first and second generation COMB potentials use an extended Lagrangian method developed by Rick et al. 49 Yu et al. 12 modified the method by adding a damp ing term for faster convergence The advantage of this Lagrangian method over solving the NxN matrix self consistently is it scales wit h O(N) instead of O(N N) On the other hand, the disadvantage of this method is it requires a fictional charge mass, a charge integration time step and an
37 arbitrary damping term, which do not correlate well to physical aspects of the material In particu lar, the charge mass in COMB potentials is 0.06 u which is two orders of magnitude larger than the actual mass of an electron. Additionally, the arbitrarily chosen convergence t olerance, usually set to 0.0001 0.01 is a trade off between charge equilibrat ion efficiency and precision, which should be chosen carefully. 2. 5 I m plementation of COMB Potentials The COMB potentials developed in this work are implemented in two MD codes, an in house HELL code and an external Large scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) 50 c ode. 2. 5 .1 Implementation to In house Code The implementation of COMB potentials in the existing in house code HELL was not a trivial task, and it certainly involved a number of researchers over the years The HELL code was originally developed and written by Drs. Dieter Wolf and Simon Phillpot. It is a two dimensional parallel MD code that includes potentials for metals, semiconductors and ionic solids such as EAM 5, 6 Tersoff 8, 9 Stillinger Weber 51 and Buckingham 7 potentials. The development and implementation of the first generation COMB potential involved a tremendous amount of coding, which includes modif yin g the HELL code and adding new potential function analysis, and fitting subroutines ; f or which Drs. Jianguo Yu and Alan McGaughey contributed significantly. The HELL implementation of COMB was first coded in Fortran 77, which used static memory allocation, the drawback of which was the inefficient usage of physical memory and frequent modification of the code related to array declaration Therefore, Bryce Devine modified the code in Fortran 90/95, which utilized dynamic memory allocation, which resulted in significantly increased efficiency. In addition to the coding and development
38 of the second generation COMB potentials, Devine also developed a serial version of the code for fitting purposes since the original parallel versi on has a lower limit on simulation size dimensions and small systems are preferred for the purpose of fitting parameters. Using the serial code to develop new functions and fitting parameters for new potentials, the computational expense is greatly reduced to the HELL implementation of COMB potentials was developing and adding necessary correction functions for specific potentials, as well as developing and implementing the multi phase fitting method. Currently both parallel and s erial versions of the COMB potentials code are still used for production and development purposes, respectively. 2. 5 .2 Implementation to Large scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) Code LAMMPS is a classical MD code that is developed and distributed by Dr. Steve Plimpton at Sandia National Laboratories. It runs on single processors or in parallel using message passing techniques with a three dim ensional spatial decomposition. LAMMPS has numerous potential functions for soft materials, such as biomolecules and polymers, and solid state materials (metals, semiconductors, ionic solids), in addition to coarse grained or mesoscopic models. It has therefore steadily gained in popularity since its release in the early 1990s. The current versi on of LAMMPS is written in C++ for the high level structure and its class organization, but functions (class methods) and potentials are written in C style codes. To allow the developed COMB potentials to become more readily available to the computational materials science/chemistry/physics community, the second generation COMB potentials that are developed in this work were implemented to LAMMPS code by the author. The
39 po tential. The efficiency of the LAMMPS implementation of COMB potentials is generally 5 10% faster than the HELL implementation (in the most efficient case for HELL ). Currently the author uses the serial HELL code for potential development and LAMMPS for pr oduction and large scale simulation works 2. 5 .3 Scaling of LAMMPS Implementation of COMB Potentials The computational expense of any simulation using any method increases with system size. Generally, the more electrons/atoms/molecules/particles/ meshes th at are explicitly treated in the system, the more CPU time is required; however the increase in computational cost does not necessarily scale linearly with system size. On the other hand, when using multiple processors computational time is efficiently red uced, but the parallel efficiency also does not necessarily increase linearly with the number of processors used. The scaling of the COMB potentials on single processor, fixed size on multiple processors and scaled size on multiple processors are discussed in following sections, respectively. All calculations are carried out on a computer cluster with 2.0 GHz Intel Xeon dual quad cores with MPICH1.2.7 message passing interface 2. 5 .3.1 Scaling with s ize on s ingle p rocessor The scaling of LAMMPS implementation of COMB potential is compared to that of LAMMPS implementation of Tersoff potential and HELL implementation of COMB potential. The potentials chosen are COMB Si/SiO 2 and Tersoff Si potentials. The efficiency is calculat ed according to: (2 4 8)
40 w he re and are the CPU time of N supercells and one supercell respectively, and N is the number of supercells An efficiency of 1.0 indicates ideal scaling. The smallest system size is a Si supercell with 4 4 4 unitcells, and the largest system size is composed of 60 of those supercells. The efficiency is plotted as a function of system size in log 10 scale, as shown in Figure 2 1. All runs are carried out with fixed and neutral charges. It is seen that both LAMMPS implementation of Tersoff and COMB potentials scale linearly with the increase in system size, however HELL implementation of COMB potential has a dec reased efficiency with increased system size. The poor scaling behavior of HELL implementation of COMB potentials is due to increased size of allocated, but not used, arr ays and less efficient utilization of link cell and neighbor list techniques. 2. 5 .3. 2 Fixed size s caling on m ultiple p rocessors The fixed size scaling on multiple processors of LAMMPS implementation of COMB potentials is carried out with SiO 2 cristobalite structure. The solution time is broken down to the time taken by inter atomic calculations, charge equilibration (QEq) calculations and their sum. All calculations consisted of 100 MD steps, and charge equilibration is conducted every 10 steps to 0.001 tolerance. The efficiency is calculated with: (2 4 9) where and are the CPU time with N processors and single processor, respectively, and NP is the number of processors used. The efficiency is plotted as a function of number of processors used, as shown in Figure 2 2. T he inter atomic
41 calculation scales linearly with increased us age of processors, indicating an ideal parallel efficiency. However, QEq shows a deviation from ideal efficiency due to increased communication among CPUs. Since QEq takes ~8 0% of the CPU time, Total time is dominated by charge equilibration. Despite the deviation from ideal parallel efficiency, a reasonable overall speedup wi th increased usage of CPUs can be expected when using LAMMPS implementation of COMB potentials. 2. 5 .3. 3 Scaled size s caling on m ultiple p rocessors The scaled size scaling on multiple processors of LAMMPS implementation of COMB potentials is carried out wi th the SiO 2 potential with a fixed 1 728 atom / CPU cristobalite structure. The system size increases with the increased use of processors so that number of atoms per processor remains constant regardless of the number of CPUs used. The largest system is r un with 64 CPUs containing 110,592 atoms. The relative computational cost in t his case is the ratio of calculation time used with multiple processors over the time used for inter atomic calculations with a single processor. The relative cost is plotted as a function of CPUs, as shown in Figure 2 3. The t otal time is the sum of QEq time and Inter atomic time plus a minor contribution from communication time. By the structure of LAMMPS, t his communication time does not include the communication time taken by the QEq calculations. It is seen in the scaled size scaling problem, inter atomic calculations remains ideal parallel efficiency, while charge equilibration dominating overall solution time and generally taking 4 5 times more time than inte r atomic calculations deviates from ideal efficiency This is also due to increased communication among increased number of CPUs during QEq calculations
42 In above scaling tests, the largest system sizes considered and the most number of CPUs used do not by any means indicate the limit of the implementation of potentials 2. 5 4 Cost of COMB Potentials Compared to other Inter atomic Potentials in LAMMPS Lastly, the relative computational cost of LAMMPS implementation of COMB potentials is compared to other implemented inter atomic potentials. As shown in Table 2 1, considered inter atomic potentials includes typical two body and three body potentials, pot entials with additional longer range dispersion and dihedral terms, as well as potentials with long range Coulombics and variable charge schemes. The LJ ratio is the computational time with single processor relative to that using Lennard Jones (LJ) potenti al, therefore larger ratio numbers indicate more expensive computational cost; the potentials are listed in the order of LJ ratio. Also shown in T able 2 1 are the system considered, number of atoms modeled timestep used and physical memory used for each p otential Data are obtained from LAMMPS official website and T able 2 1 is reproduced with permission. It is seen that when three body or many body interactions are considered, the relative cost to LJ is increased by a factor of 4, including long range disp ersion or Coulombic increases the cost by approximately an order of magnitude, and further including dynamic charge equilibration increases the cost by two orders of magnitude. Charge equilibration of COMB potentials took ~230 of the computational time, while the inter atomic calculations, including long range Coulombics, took only ~50 computational time. Therefore provided the efficiency of charge equilibration is greatly increased, the overall computational cost of COMB potentials can be significantly reduced. In addition, physical memory of the computing node places a hard limit on the
43 maximum atoms allowed per processor. It is seen from T able 2 1 that the COMB potentials use far less physical memory compared to potentials with roughly the same order of LJ ratio, which indicates that compared to ReaxFF and eFF potentials larger system sizes are allowed when using COMB potentials 2.6 Summary The COMB potentials developed in this work equilibrate atomic charges based on the electronegativity equalization principle, and the charge of a particular atom is equilibrated according to its bonding environment. These developed COMB potentials are implemented in HELL and LAMMPS 50 codes, and the efficiency and scalability of the LAMMPS implementation is shown to be better than that of the HELL code. This is mainly because HELL code is designed for parameterization and development of new potentials and is not optimized for computati onal efficiency.
44 Table 2 1. Relative computational cost of selected inter atomic potentials available in LAMMPS code. Re produced with permission from Steve J. Plimpton, http://lammps.sandia.gov Potential System Atoms Timestep Memory LJ ratio Lennard Jones LJ liquid 32,000 0.005 tau 12 Mb 1.0 Embedded Atom Method bulk Cu 32,000 5 fs 13 Mb 2.4 Tersoff bulk Si 32,000 1 fs 9.2 Mb 4.1 Stillinger Weber bulk Si 32,000 1fs 11 Mb 4.1 Embedded Ion Method NaCl crys. 32,000 0.5 fs 14 Mb 6.5 CHARMM+PPPM solvated protein 32,000 2 fs 124 Mb 13.6 Modified EAM bulk Ni 32,000 5 fs 54 Mb 15.6 AIREBO polyethylene 32,640 0.5 fs 101 Mb 54.7 ReaxFF/C PETN crys. 32,480 0.1 fs 976 Mb 185 COMB SiO 2 cyrs. 32,400 0.2 fs 85 Mb 284 eFF H plasma 32,000 0.001 fs 365 Mb 306 ReaxFF PETN crys. 16,240 0.1 fs 425 Mb 337
45 Figure 2 1. Scaling with size on single processor for LAMMPS Tersoff, LAMMPS COMB and HELL COMB. Computation efficiency plotted as a function of system sizes. 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1 10 100 Efficiency Number of supercells (log scale) Scaling on single processor LAMMPS_Tersoff LAMMPS_COMB HELL_COMB
46 Figure 2 2. Fixed size scaling of COMB potential. Computation efficiency plotted as a function of number of CPUs used for a fixed system size problem. 0.0 0.2 0.4 0.6 0.8 1.0 1.2 0 5 10 15 20 25 30 35 Efficiency Number of CPUs Fixed size scaling of LAMMPS COMB Total QEq Interatomic
47 Figure 2 3. Scaled size scaling of COMB potential. Computation efficiency plotted as a function of number of CPUs used for a fixed size per CPU problem. 0 1 2 3 4 5 6 7 0 10 20 30 40 50 60 70 Relative cost Number of CPUs Scaled size scaling of LAMMPS COMB Total Qeq Interatomic
48 CHAPTER 3 DEVELOPMENT OF COMB SI/SIO 2 POTENTIAL AND ITS AP PLICATION S TO AMORPHOUS SILICA AND SI/SIO 2 INTERFACES Silicon is ubiquitous in microelectroni c devices and ultra large integrated circuits (ULSIs); silicon dioxide (SiO 2 silica), in particular the amorphous phase was used as the gate oxide for decades. The most stable crystalline phase at ambient conditions is quartz (space group P3 2 21 ). There are three other ambient pressure, high quartz ( P6 2 22 tridymite ( P6 3 /mmc cristobalite ( Fd m ). cristobalite melts at 2023K. Silica crystallizes into two crystalline polymorphs at higher pressures: coesite ( C2/c ) above 2 GPa and stishovite ( P4 2 /mnm ) above 9 GPa. 52, 53 The cristobalite ( P4 1 2 1 2 ) phase is me tastable relative to quartz but is lower in energy than all high temper ature and high pressure phases. Recently Yu et al. developed a COMB potential that allows dynamic charge equilibration for Si/SiO 2 systems 12 (referred to a s COMB07 hereafter). While COMB07 does a good job of describing the phase order among silica polymorphs and includes a variable charge scheme the elastic constants of various silica phases were poor, and amorphous silica was not well described. In this wo rk, therefore, a second generation form of this potential (referred to as COMB10 hereafter) that includes a modified functional form alism for the potential and a new parameterization is developed The details of COMB10 potential are described in Section 2.2. COMB10 is applied to both Si/SiO 2 interfacial systems and amorphous silica; it not only describes the phase order, structural properties, and mechanical properties of the various silica polymorphs, but it also does a good job of describing amorphous s ilica, atomic oxygen adsorption
49 energies, and Si/SiO 2 interfaces. The developed COMB10 Si/SiO 2 potential is also applied to investigate the properties of Si nanocrystals embedded in amorphous silica and to carry out large scale nanoindentation simulations of Si/SiO 2 systems, including quartz and amorphous silica. This work is only made possible with the contribution from and helpful discussions with numerous researchers. Bryce Devine not only assisted tremendously in the development and implementation of the parameter fitting and property analysis subroutines, but also determined the oxygen parameters that describe O O interactions. Jeffery Hawkins and Aravind Asthagiri contributed in the first principles calculations of atomic oxygen adsorption energies, and Kai Nordlund and Flyura Djurabekova supplied the Si nanocrystals embedded in amorphous silica structures. Scott Warren assisted in writing a program that generates various geometries of indenters used in nanoindentation simulations. In addition, t he au thor also gratefully acknowledge s the contributions of Kai Nordlund and Flyura Djurabekova for supplying the amorphous silica structures used in finalizing the parameter set. 3. 1 Potential Development and Properties of SiO 2 Polymorphs For the COMB10 Si/SiO 2 potential, t he short range parameters for Si are the same as the Tersoff potential for Si 9 The author, with assistance from Devine determined the O O and Si O interaction parameters; t he parameterization is given in Table A 1 Table 3 1 compares the properties of cristobalite phases calculated with the COMB10 potential with experimental data and DFT (with Perdew Burke Ernzerhof (PBE) exchange correlation functional 54 ) results, COMB07 predictions, and values obtained using two other empirical potentials, TTAM 55 and BKS 56 The properties given cristobalite are predicted. The table
50 quartz reasonably well compared to the TTAM and BKS potentials, while showing substantial cristob alite show deviations from the experimental values, COMB10 predicts properties with better materials fidelity than COMB07. cristobalite given by COMB07 resulted from the strong coefficients 4.5 eV and 9.5 eV respectively, on the Si O Si and O Si O bond bending terms ( Section 2.2. 4.1, Equation s 2 3 9, 2 4 0) and from the strong (20 eV) additional short range repulsion ( Section 2.2. 4.3, Equation 2 4 3) Strong penalization for the deviation of the above bonds from the ideal bond quartz) was necessary for COMB07 to reproduce the correct phase order for silica polymorphs, though at the expense of elastic constants. The bond bending terms in COMB10 are significan tly weaker 2.6 eV and 0.31 eV respectively, than in COMB07. Additionally, describing the Coulombic interactions with integrals over spherical charge densities provided enough repulsion between charged atoms (Section 2.2. 2, Equation s 2 3 5, 2 3 6); as a resu lt no additional short range repulsion is needed for COMB10 In turn, this leads to elastic constants that are much closer to experimental values than COMB07 It should be noted that other cristobalite better than those given in the COMB10 column in Table 3 1 ; however, they were inferior for other polymorphs, in some cases to the point of giving an incorrect phase order. Table 3 2 gives the phase order of silica polymorphs by quantifying th e cohesive quartz (in units of eV/SiO 2
51 unit) as calculated with COMB10 and COMB07, and as obtained from experiments, first principles calculations, and from other empirical potentials. The table sho ws that COMB10 gives the correct phase order for two high quartz and quartz. However, the third high tr idymite by COMB10, which contradicts DFT literature values. Nevertheless, COMB10 correctly quartz, and the right phase order for the high pressure phases. The significant overestima tion in the cohesive energy of stishovite results from the bond bending terms in the potential, since stishovite has 6 fold coordinated Si atoms and the Si O Si and O Si O bond angles are substantially different from all the other polymorphs considered. Th e finite temperature properties are also of great interest to materials scientists. The COMB10 potential predicts the thermal expansion coefficients of quartz in the a and c directions to be 8.810 6 K 1 and 15.010 6 K 1 respectively; both are within + 13% of the experimental values. 57 quartz from COMB10 is approximately larger than the experimental value by a factor of two. 3. 2 Propert ies of Amorphous Silica and its Defect Formation Energies A physically appropriate description of amorphous silica is also extremely important for the atomistic modeling of microelectronic devices and ULSIs. COMB07 was not able to describe amorphous silica the obtained high temperature amorphous silica melt re quartz upon cooling due to strong bond angle penalization In this work, we use the developed COMB10 potential for Si/SiO 2 to obtain amorphous silica and investigate its defect for mation energies.
52 3. 2 .1 Structural Properties of Amorphous Silica To obtain amorphous silica cristobalite is heated from 300 K up to 3000 K for 100 ps with a constant temperature and constant pressure ensemble, allowing the simulation box size and shape to change. The silica melt is then maintained at 3000 K at a constant box volume for 100 ps to ensure complete melting. The silica melt is then quenched (using a damped force minimization method) followed by annealing at 300 K for 50 ps allowing the cell size and symmetry to relax. The final density and cohesive energy of the amorphous silica thus obtained are 2.458 g/cm 3 and 20.50 eV/SiO 2 ; the average c harge on the Si and O atoms are +2.898 and 1.449, respectively. The averaged values of the Si O bond length, O O distance in a particular SiO 4 tetrahedral unit, and O Si O bond angle thus obtained are 1.63 2.71 and 107.2, respectively. These are in fair agreement with the experimental 58 values of 1.6 2.6 and 109.28. COMB10 predicts a very wide distribution of Si O Si bond angles (110 160), with an average value of 135.4; experiments yield angles in the range 130 160, with an average value of 144 59 The radial distribution function and bond angle distribution of the amorphous silica is shown i n Figure 3 2. 3. 2 .2 Point Defect Formation Energies of Amorphous Silica Table 3 3 compares the defect energies of the cation and anion vacancies and interstitials of amorphous silica calculated with COMB10 to DFT LDA results 60, 61 The defect formation energies are calculated in the same way as in Shan [ 13 ] These defect formation energies are not part of the fitting set of physical properties, but are predictions. The table indicates th at the defect formation energies for a Si interstitial and an O vacancy from the COMB10 potential are slightly larger than those from the DFT calculations, while the formation energies of a Si vacancy and an O interstitial are
53 significantly overestimated. Nevertheless, the general trend of the defect formation energies is captured, and the oxygen interstitial is predicted to have the smallest formation energy, which is consistent with the results of the DFT calculations. Also shown are the estimated defect formation energies of Schottky defects, and cation and anion Frenkel pairs. These results indicate that the energy required to create a Schottky defect is approximately the same as that of one Si and two O vacancies added together. Likewise, the formation energies of the cation and anion Frenkel pairs are very close to the sum of their constituent point defects, indicating negligible association energies of vacancy and interstitial point defects. 3. 3 Atomic Oxygen Adsorption Energies on Reconstructed Si (0 01) Surface The adsorption of oxygen atoms and molecules on Si surfaces plays a crucial role to the Si oxidation processes. To evaluate the ability of the COMB10 potential to describe the energetics of oxygen adsorption, we characterize the adsorption ener gies of atomic oxygen on a 2 1 reconstructed Si (001) surface with DFT calculations and MD simulations using COMB10 potential. The DFT calculations are performed using the Vienna ab initio simulation package (VASP). 62 We use the projector augmented wave (PAW) pseudopotentials 63, 64 provided in the VASP database. Calculations have been done using the generalized gradient approximation 65 Perdew Burke Ernzerhof (GGA PBE) exchange correlation functiona l 54 including spin polarization. A plane wave expansion with a cutoff of 400 eV is used, and the total energy calculations are done using a block Davidson iteration method for electronic relaxations, accelerated using Methfessel Paxton 66 Fermi level smearing with a Gaussian width of 0.2 eV. The positions of the atoms are relaxed until the forces on all unconstrained atoms are less than 0.03 eV/. All calculations performed on an eight layer slab with a 22 surface
54 unit cell. The adsorbed O atoms and the top five Si layers are permitted to fully relax while the lower three Si layers are held fixed. A vacuum region of ~12 ensures that the slab does not interact with its periodic image in the surface normal direction. We use a 661 Monkhorst Pack k point mesh 67 We have tested larger plane wave cutoffs and slab thickness to confirm convergence of the DFT adsorption energies. The lattice parameter for Si is found to be 5.469 using DFT versus an experimental value of 5.43 The DFT calculations are done wit h spin polarization and the minima can take on singlet or triplet configurations and we have explored both spin states for all the adsorption sites. The COMB potential formalism is not able to produce information on the spin states of the minima and we foc us only on the adsorption energies. The DFT results are in good agreement to earlier B3LYP DFT cluster studies of O adsorption on Si(001). 68, 69 The MD s imulations are carried out with the COMB10 potential as implemented in the LAMMPS 50 program. The Si supercell, whose size is ~ 30 30 15 with a vacuum space of 25 in the  direction, consists of 768 atoms and two free surfaces. 16 O atoms are adsorbed on each Si surface that results a 0.5 monolayer surface coverage. For purposes of comparison, the surfaces are quenched (force and energy minimized) and then equilibrated with charge equilibration at a constant temperature of 0.1K for at least 1ps until minimal fluctuation in total energy is observed. The adsorption energies, E ads are calculated in the standard manner: E ads = (E Si + E O ) E Si O where E Si E O and E Si O are the energies of the isolated Si surface, O atoms, and Si su rface with O atoms adsorbed, respectively; th e results are shown in Table 3 4 where we have reported the most favored O configuration at the atop, bridge,
55 and backbond sites. It is seen from DFT PBE calculations that the atop and backbond adsorption sites have the smallest and largest adsorption energies, respectively, indicating that the atop site is the least favorable adsorption site on the 2 1 reconstructed Si (001) surface, while the backbond site is the most favorable; the bridge site has an interm ediate adsorption energy. Adsorption energies from MD simulations using COMB10 potential are comparable to DFT calculations, with a difference of less than 6% even in the largest case; these differences are shown in parentheses. Figure 3 3 shows the relaxed structures of oxygen atoms adsorbed at the atop site, the bridge site and the backbond site, from left to right respectively, on the Si (001) surface from (a) DFT PBE calculations and (b) the COMB10 potential. Atomic charges from DFT cal culations are Bader 70 charges. Although the charges of the adsorbed O from COMB10 potential are only ~65 % of that from DFT calculations, the COM B10 potential is capable of predicting the change in the distribution of atomic charges when an O atom is adsorbed on the Si surface in a quite reasonable manner. The average Si O bond length between the adsorbed O and the neighboring Si atoms for the atop bridge and backbond sites obtained from DFT PBE calculations are 1.337 1.681 and 1.535 respectively, compared to 1.478 1.727 and 1.665 obtained from MD simulations. It is seen that not only the adsorption energies from COMB10 correspond w ell with DFT calculations, but also the predicted structures are reasonably comparable. A systematic investigation of the adsorption energies of molecular oxygen on the Si surfaces is reserved for a future study.
56 3. 4 Si/SiO 2 Interfaces 3. 4 .1 Si/ Q uartz I n terface The performance of silicon based electronic devices such as MOS/CMOS field effect transistors is known to be significantly affected by the Si/SiO 2 interface. Computational studies of such interface mostly rely on DFT calculations or MD simulations with empirical potentials such as BKS 56 TTAM 55 and ReaxFF 11 While accurate, DFT calculations are computationally expensive and are limited t o systems with less than several hundred to a few thousand atoms; it is therefore difficult to exclude unrealistic mechanical forces and/or polarization. Compared to fixed charge empirical potentials (BKS, TTAM), variable charge potentials allow the adjust ment of atomic charges in a response to changing atomic environment, such as across an interface. To characterize the ab ility of COMB10 potential to model a Si/SiO 2 interface, we quartz(010) interface. The inter facial system quartz slab both made of 10 510 unstrained conventional unit cells. To make the interface fully coherent, initial 8.3% and 2.3% strains are applied to the  and  directions, respectively, to the Si slab. To find the most stable configuration of the interface with the lowest energy, the O terminated quartz slab is shifted by an integer number (0, 1, 2, 3) of the lattice parameters with respect to the Si slab in both  and  directions, resulting in a total of 16 different interfacial configurations. Force and energy minimization is carried out on these configurations while simultaneously allowing the atoms to move, charges to vary, and the simulation box size to change. The interfacial configuration with the lowest energy is shown in Figure 3 4 (a) and (b) as snapshots of the interfacial structure after
57 force/energy minimization viewing along  and  directions, respectively, and color coded by t he coordination number (black = 2 fold, white = 4 fold). The interface stays coherent and contains Si coordination defects half of the interfacial Si atoms are two fold coordinated and the other half is four fold coordinated, while all O atoms are two fo ld coordinated. Comparing to other minimized configurations, this has the least number of coordination defects, thus making the quartz slabs are bonded through Si O bonds at the interface, with every two interfacial O atoms bonded to one interfacial Si atom. The interfacial Si O bond length averaged to be 1.89 which is ~15 % larger quartz. The minimized box sizes in  and  directions correspond to a slightly contracted S quartz slab compared to the equilibrium lattice parameters. The Si slab is contracted 6.8% and 1.5% in  and  quartz slab is expanded +1.6% and +0.8% in  and  directions. The min imized interface is then heated at constant pressure, allowing the box size to change, with the temperature gradually brought up to 500K in 12 ps, followed by another 8 ps of equilibration at 300K. The snapshots of the final equilibrated interfacial system at 300K are shown in Figure 3 5 (a) and (b), viewing along the  and  directions, respectively, and are color coded by the atomic charges with charge values indicated by the color bar (white is positive, gray is neutral and black is negative). Th e equilibrated interfacial system has an average Si O bond length of 1.80 and the number of Si O bonds also increases, with every three O atoms at the interface bonded to two interfacial Si atoms both indicating a stronger bonding between the two slab s.
58 The coordination defects at the interface after equilibration have changed from under coordinated Si atoms to over coordinated O defects. The interfacial Si atoms are all four fold coordinated, while two of the three interfacial O atoms are three fold c oordinated; these over coordinated O atoms have two shorter than average and one longer than average Si O bonds. The distribution of charges of Si and O atoms across quartz are illustrated in Figure 3 6 The atoms further away from the interface i quartz and Si slabs have normal bulk atomic charges, while the interfacial Si and O atoms have intermediate and more negative charges, respectively. The charge distribution at the interface indicates the donation of electrons from Si atoms to O atoms f orming strong Si O bonds. 3. 4 .2 Si N anocrystals E mbedded in A morphous S ilica Si/SiO 2 systems and interfaces are widely used in modern optoelectronic industries, particularly systems such as Si quantum dots (nanocrystals) embedded in amorphous silica matrix. These structures are formed by ion beam implantation of Si ions into amorphous silica separating the SiO 2 phase into small Si nanocrystals with diameters dependent upon post implantation annealing temperatures. 71 The optical properties and the nature of the optical responses of these systems are believed to arise either from the splitting of the energy levels of the Si nanocrystals due to the quantum confinement effect 72, 73 or from the radiative recombination centers for excitons at the Si nanocrystals and amorphous silica interfaces 74 76 in particular double Si=O ( silanone ) bonds 75 In this work, therefore, we apply the developed COMB10 Si/SiO 2 potential to examine the structure of Si nanocrystals embedded in an amorphous silica matrix (SiNC /a SiO 2 ). The potential is not able to reproduce quantum confinement effects;
59 nevertheless it is able to investigate the interfacial properties and the charge distribution of the system. The SiNC/a SiO 2 structures were provided by Flyura Djurabekova and Ka i Nordlund. 77 quartz crystal structure wi th high energy primary knock on atoms method to create amorphous silica matrix followed by carving a spherical hole out of the matrix and replacing the void with crystalline Si spheres. The system is equilibrated with the COMB10 potential, and a snapshot after equilibration is given in Figure 3 7. The simulations predict that interfacial Si atoms have charges around +1.5e, which is intermediate between neutral charged Si atoms in the Si crystal and +2.9e charged Si atoms in a SiO 2 Silanone bonds can be fo und at the Si/SiO 2 interfaces, shown in Figure 3 8, identified by a shorter Si O bond length and larger fractional charges on the Si and O atoms. Initial results indicate the ability of COMB10 potential to be applied to such SiNC/a SiO 2 systems and furthe r in depth investigation is required to draw conclusions of the properties of Si nanocrystals embedded in amorphous silica matrixes. 3. 5 Nanoindentation of Si/SiO 2 Systems Indentation test ing is one of the most commonly applied methods of probing the mecha nical properties of materials. 78 84 In such a test, a hard tip, also known as an indenter, with known mechanical properties, is pressed with an applied load into a sample whose mechanical properties are unknown. The load is applied, sometimes held, and then removed. The area of the residual indentation in the sample is measured t o evaluate the hardness, H by dividing the maximum load, P by the residual indentation area, A r or (3 1)
60 while the slope of the initial portion of the unloading curve provides the stiffness S : (3 2) where h is the depth of the indenter at the onset of unloading. The stiffness can be used to calculate the reduced modulus E r of the substrate with the following equation : (3 3) where A is the contact area between the indenter and the substrate and can be approximated with the following equation for a spherical indenter 85 : (3 4) where is the change in heig ht of the surface before and after indentation. Nanoindentation, which is indentation testing in the submicrometer or nanometer to 400 nm, respectively. This hardness te sting method is an important tool and is considered a standard, for both scientific research and m aterials characterization Since nanoindentation tests measure the mechanical properties of volumes so small that they can be made defect free, it probes the fundamental processes that initiate deformation. This is substantial for understanding the ideal strength of materials, which is the strength that ideal, defect free crystals can sus tain. Although powerful microscopy techniques such as friction force micro scopy (FFM), atomic force microscopy (AFM) and in situ transmission electron microscopy (TEM) are capable of providing detailed images with atomic scale resolution, it is still challenging to interpret nanoindentation experiments at a fundamental level bec ause of the very complex stress profile and atomic movements generated in the sample near the
61 indenter. Analytic models and atomistic computational simulations have been very helpful in playing the role of characterizing and understanding the processes und erlying the nanoindentation response. 86 92 The mechanical properties of materials being used in microelectronic devices, such as oxides and carbides, at the submicron range are still not yet fully understood. Nanoindentation MD simulations have been applied to investigate the mechanical properties of these materials including SiC, 93 98 ZnO 99 Al 2 O 3 100 and MgO 101 as well as that of interfac es. 102 Si/SiO 2 systems are widely used in electronic devices however, t o currently no published results on MD simulations of nanoindentation of Si/SiO 2 systems/interfaces can be found in the literature This is most likely due to the limitation of an appropriate empirical potential function. Therefore, here we apply the developed COMB10 Si/SiO 2 potentials to study the response of the Si/SiO 2 systems/interfaces during nanoindentation simulations. The Si/SiO 2 systems include Si/a SiO 2 quartz interfaces, as shown in Fig s 3 9 (a) and (b), respectively. The Si/SiO 2 systems consist of Si and SiO 2 (a SiO 2 quartz) layers with thicknesses of 27 and 1 2 respectively, and the dimensions normal to the indentation direction are around 50 The interfacial systems are pressed with a rigid hemi spherical Si indenter with a 33 diameter of curvature at a constant displacement rate of 1m/s. Nanoindentation simulatio ns are carried out at 300 K with Nose Hoover thermostat 20 22 applied to bottom half of the Si slab while the rest of the atoms are allowed to evolve without any additional constraints. During the constant di splacement loading of the indenter into the Si/SiO 2 systems, the total load reflected to the indenter is recorded as a function of the displacement of
62 the indenter. The load displacement curves are shown in Figs. 3 10 (a) and (b) for Si/a SiO 2 quartz interfaces, respectively, and it is seen that the loading stage resulted in plastic deformation to the Si/SiO 2 systems. The hardness and reduced modulus of the Si/SiO 2 systems from this nanoindentation tests are given in Table 3 5. The Si/ quartz interfacial system has a larger hardness and reduced modulus than that of the Si/a SiO 2 system, which is consistent with experimental findings conducted to quartz and a SiO 2 substrates. 80 3. 6 Summary This new COMB10 poten tial captures many of the physical properties of silica polymorphs and amorphous silica; it can also be applied to model Si/SiO 2 interfacial systems including the properties of Si quantum dots embedded in amorphous silica and mechanical testing by means o f nanoindentation of Si/SiO 2 systems It should be noted that TTAM and BKS potentials do an overall better job than COMB10 in describing the properties of several SiO 2 polymorphs. However, they are unable to consider either pure silicon or Si/SiO 2 interfac es Thus, the COMB10 potential for Si/SiO 2 should prove to be a useful new tool for the computational toolbox and an effective method for carrying out large scale atomistic simulations of systems of technological importance.
63 Table 3 1 Comparison of the cristobalite calculated by the COMB10 potential with values obtained from experiments, first principles calculations (in parentheses), COMB07 and two other empirical potentials (TTAM and BKS) quartz Properties Expt. (DFT) COMB10 COMB07 12 TTAM 55 BKS 56 a 0 () 4.916 4.856 4.978 5.02 4.94 c 0 () 5.405 5.316 5.307 5.55 5.45 E c (eV/SiO 2 ) 19.23 ( 22.81) 20.63 19.40 17.85 B (GPa) 39 64 134 33 41 G (Gpa) 47 39 116 36 45 C 11 (GPa) 87 99 238 71 91 C 12 (GPa) 7 5 61 9 8 C 13 (GPa) 13 38 91 12 15 C 14 (GPa) 18 10 0.7 14 18 C 33 (GPa) 107 111 242 91 107 C 44 (GPa) 58 42 167 40 50 q Si (e) (4.0) 2.92 1.70 2.4 2.4 cristobalite Properties Expt. COMB10 COMB07 12 TTAM 103 BKS 103 a 0 () 4.98 5.02 4.94 c 0 () 6.94 6.94 6.71 E c (eV/SiO 2 ) 19.20 20.57 19.15 C 11 (GPa) 59 137 454 48 65 C 12 (GPa) 4 18 27 6 7 C 13 (GPa) 4 43 62 4 1 C 33 (GPa) 42 118 396 35 38 C 44 (GPa) 67 55 186 58 70 C 66 (GPa) 26 29 149 20 28 q Si (e) 2.89 1.75
64 Table 3 2 quartz (in units of eV/SiO 2 unit) and of six polymorphs quartz calculated with the COMB10 potential. The values are compared to published results from experiments, first principles calculations and other empirical potentials. (eV/atom) Exp COMB 10 DFT 104 COMB 07 12 TTAM 55 ReaxFF 11 Yasukawa 25 Quartz 19.23 20.63 25.964 19.33 22.2 19.275 Quartz +0.051 +0.108 +0.026 +0.237 +0.073 0.03 Tridymite +0.635 +0.034 +0.615 +0.258 0.018 +0.00 Cristobalite +0.054 +0.500 +0.060 +0.582 +0.213 0.765 Cristobalite +0.03 +0.049 +0.025 +0.036 +0.186 +0.003 2.865 Coesite +0.03 +0.082 +0.012 +0.267 +0.018 Stishovite +1.196 +0.104 +0.837 Table 3 3 Defect formation energies in amorphous silica predicted with the COMB10 potential compared with those from published DFT calculations Defect DFT (eV) 60 DFT (eV) 61 COMB10 (eV) Si vacancy 5.69 4.78 16.51 Si interstitial 12.84 14.07 18.05 O vacancy 6.20 5.71 7.07 O interstitial 1.86 1.56 5.44 Schottky --31.54 Cation Frenkel --34.02 Anion Frenkel --11.25 Table 3 4 Atomic oxygen adsorption energies (eV) at three different adsorption sites on a 2 1 reconstructed Si (001) surface from DFT PBE calculations and COMB10 potential. Atop Bridge Backbond DFT PBE 5.42 6.35 6.48 COMB10 5.71 (+5.3%) 6.10 ( 3.9%) 6.24 ( 3.7%) Table 3 5. Hardness and reduced modulus of Si/ quartz and Si/a SiO 2 systems obtained from nanoindentation simulations using COMB Si/SiO 2 potential developed in this work compared to experimental values. Type of SiO 2 layer H (GPa) E r (GPa) Quartz 45.1 (~105 a ) 221.0 (125 a ) Amorphous silica 34.0 (~80 a ) 202.6 (71 a ) a Ref. [ 80 ]
65 Figure 3 1. The cohesive energy, in units of eV/SiO 2 unit, as a function of unit volume curves for SiO 2 polymorphs from COMB Si/SiO 2 potential developed in this work.
66 A B Figure 3 2. A morphous silica obtained using COMB10 Si/SiO 2 potentia l developed in this work A ) Radial distribution function. B ) B ond angle distribution Separation, r () O Si O and Si O Si bond angles ()
67 A B Figure 3 3 The final relaxed structures of oxygen atoms adsorbed at the atop site, bridge site and backbond site, from left to right respectively, on a 2 1 reconstructed Si (001) surface from A ) DFT PBE calculations, and B ) MD using the COMB10 potential. Larger atoms are Si and the smaller ones are O; atoms are gray scale color coded by the charges with the c olor bar indicating the charge values.
68 A B Figure 3 4 quartz interfacial structure after force and energy minimization, viewing along A )  and B )  directions. Larger spheres are Si atoms, smaller ones are O atoms. The atoms are gray scale color coded by the coordination number, with black and white being two and four fold coordinated, respectively.
69 A B Figure 3 5 quartz interf acial structure after equilibrating at 300K for 20ps, viewing along A )  and B )  directions. Larger spheres are Si atoms, smaller ones are O atoms. The atoms are gray scale color coded by the atomic charges, with the charge values indicated by th e color bar.
70 Figure 3 6 quartz interface. Length in Y direction ( )
71 Figure 3 7 Snapshot of a Si nanocrystal, 6 nm in diameter, embedded in a SiO 2 matrix from MD simulations using COMB Si/SiO 2 potential. Atoms are color coded by the charges indicated by the color bar (red=positive, cyan=neutral, blue=negative)
72 Figure 3 8 Silanone bonds found at the Si/SiO 2 interface of the Si nanocrystal embedded in a SiO 2 matrix. Atoms are color coded by the charges indicated by the same color bar as Figure 3 6.
73 A B Figure 3 9 Snapshots of A ) Si/a SiO 2 and B ) quartz systems from nanoindentation simulations using COMB Si/SiO 2 po tential developed in this work. Hemispherical rigid indenter on the top, SiO 2 layer in the middle, Si slab on the bottom. Atoms are color coded by the charges with the same color bar
74 A B Figure 3 10 Load displacement curves of A ) Si/a SiO 2 and B quartz systems from nanoindentation simulations. Loading and unloading rates are 1 m/s; loading curves shown in red, unloading curves shown in blue.
75 CHAPTER 4 DEVELOPMENT OF COMB CU/CU 2 O POTENTIAL AND ITS AP PLICATION TO CU/SIO 2 INTERFACES Interfaces between metals and oxides are important because of their ubiquitous presence in numerous materials structures, including microelectronic devices. Copper, one of the most used inter connect material in silicon based microelectronics and related devices, has low electrical resistivity and high electromigration resistance. 105, 106 It has been reported that high purity Cu films bind well to the silica substrat es if there are no hydroxyl groups at the Cu/SiO 2 interface. 107, 108 A ma jor issue, however, is the interaction between Cu and SiO 2 at the interface; this interaction can lead to the formation of oxidized Cu leading to the diffusion of Cu ions through the SiO 2 layer, which results in the degradation of the dielectric layer. 109, 110 Computational studies of Cu/SiO 2 interfaces using first principles density functional theory (DFT) calculations have provided valuable insights into the structure and energetics of these interfaces. For example, Nagao et al. 108 characterized the structural, electronic, and adhesive properties of interfaces between fcc cristobalite (001) with different types of interface terminations 3, 4 and predicted that adhesion at the Cu/SiO 2 i nterface is the strongest in the most oxygen rich case. The calculated adhesion energy was consistent with the value obtained experimentally by Kriese et al. 111 through nanoindentation experiments and Pang et al. 112 through delamination experiments. The work reported in this chapter utilizes classical atomistic simulations with combined COMB potentials for Si/SiO 2 16 and Cu/Cu 2 O 17 cristobalite (001) Q) and Cu (001)/amorphous silica (referred as Cu/a SiO 2 ) interfaces. The potentials for Si/SiO 2
76 and Cu/Cu 2 O follow the COMB10 formalism, discussed fully in Section 2.2. These simulations allo w us to extract the structural and adhesive properties, and charge transfer across the interface. The fidelity of the MD simulations and COMB potentials are validated against the results of DFT calculations or experimental values for the cases where they a re available. 4. 1 Potential Development and Properties of Cu and Cu 2 O The parameterization of the Cu 2 O and Si Cu interactions considered both the C interfaces and Cu 2 O bulk properties. The Si Cu mixing scaling factors shown in Ta ble A 2, are only applied to Si Cu bonds when the Si atom is also bonded to an O atom; the scaling factors are set to 1.0 for all other Si Cu bonds. However, Cu O scaling factors are applied to Cu O bonds regardless of the additional bonding of the Cu atom The properties of Cu 2 O from the COMB potential compared to experimental and first principles calculations are given in Table 4 1 It is seen the lattice parameter of Cu 2 O with this potential shows a ~3% deviation from first principles calculations but the cohesive energy is well reproduced. Some other sets of parameters and mixing rule scaling factors actually give better Cu 2 O properties than the set given in Table A 1 ; however these parameters and scaling factors give worse adhesion energy C i nterfaces. The properties of Cu and SiO 2 from the COMB potentials can be found in Yu [ 40 ] and Shan [ 16 ] (also Chapter 3) respective ly. 4. 2 Adhesion of Cu/SiO 2 interfaces 4. 2 .1 cristobalite Interfaces To provide reference results against which to compare the results of our atomistic simulations, we perform first principles density functional theory calculations of the C interfaces within the generalized gradient approximation (GGA), 65 using the
77 Perdew Burke Ernzerhof (PBE) exchange correlation functional, 54 as implemente d within the Vienna ab initio simulation package (VASP). 62 We use plane wave basis sets with a 500 eV energy cutoff, and projector augmented wave (PAW) pseudopotentials 63, 64 for Si, O, Cu and H. The convergence criteria are set at 1 .010 4 1 for energies and forces, respectively. We also make use of a Fermi distribution 66 smearing with a temperature of k B T ~ 0.2 eV. For purposes of comparison, we also carry out local de nsity approximation (LDA) calculations with a 394 eV energy cutoff, and used ultrasoft pseudopotentials for Si, Cu and H, and PAW potentials for O. The rest of the computational setups are kept the same as in the GGA calculations. To enable us to make a co mparison with the work of Nagao et al. 108 we examine the same fcc cristobalite (001) interface that they considered in their work. cristobalite slab consist s of nine SiO 2 layers in the  direction with a 22 surface unit cell. The dangling O atoms at the bottom are hydrogen terminated and are fixed at their bulk positions to reduce size effects, while all other atoms are allowed to fully relax. The Cu sla b consists of six (001) layers with a 22 surface unit cell, whose orientation is rotated 45about the  axis so that the C interface is varied to mimic oxygen rich or oxygen C interfaces: OO terminated (Cu/ C:OO), O terminated (Cu/ C:O) and Si terminated (Cu/ C:Si) C inter facial system is ~27 The Monkhorst Pack k point mesh 67 is 441 for the 9.92 9.92 45 supercell (with an 18 vacuum region to prevent interacting with its own periodic image ). The relaxed
78 structures of the Cu/ C:OO, Cu/ C:O, and Cu/ C:Si interfaces with GGA PBE viewed along the  and  directions are shown in Fig s 4 1 (a) and (b), respectively, which are quite consistent with that obtained with LDA shown in Nagao [ 108 ]. In particular, the relaxed Cu O and Cu Si bond length at the interface from our GGA PBE calculations are ~1.93 and ~2.47 respectively, compared to ~1.90 and ~2.40 respectively from LDA. 108 These differences in Cu O and Cu Si bond lengths are consistent with the general notion that LDA calculations tend to overestimate bonding energy compared to GGA PBE calculations. C interfaces for our atomistic simulations using COMB potentials are those used in the DFT calculations, except that the simulation cells have a larger area, and are thicker. Force and energy minimization are applied to find the lowest energy interfacial configurations for the Cu/ C:OO Cu/ C:O and Cu/ C:Si interfaces, as illustrated in Figs. 4 2a and 4 C interfaces from COMB are strikingly similar to those obtained from our GGA PBE calculations. The minimized Cu O and Cu Si bond length at the interfaces from ou r atomistic simulations are ~1.86 and ~2.32 respectively, which are ~5 % smaller than those from GGA PBE calculations. The smaller Cu O and Cu Si bond lengths from COMB potentials compared to DFT calculations come from generally smaller fitted lattice parameters for bulk Cu 2 O phase as presented in Table 4 1 and Shan [ 37 ] and SiO 2 polymorphs as presented in Table 3 1 and in S han [ 16 ] The adhesion energy, W is calculated from: (4 1) where and are the energies of the isolated SiO 2 and Cu slabs and the C interfaces, respectively, and A is the surface/interface area. The calculated
79 adhesio n energies are given in Table 4 2 The COMB potentials reproduce the same trends for the adhesion energies as the DFT calculations: oxygen rich interface exhibits the strongest adhesion, followed by intermediate and oxygen lean interfaces. These are consistent with the notion that Cu films bind well to oxygen rich SiO 2 substrates. 107, 108 The reason of which is more Cu O bonds formed at the interface due to larger numbers of oxidized Cu atoms. The fidelity of the combined COMB potentials is thus validated against the results of DFT calculations, and it is shown that COMB potentials are adequate to model Cu/ C interfaces. Therefore in the following results and discussion we apply COMB potentials to investigate more Cu/SiO 2 interfaces. The intermediate oxygen concentration interface, that is Cu/ C:O, c onsidered above is constructed by removing one of the four O atoms from the surface Si from the Cu/ C:OO interface. In this case every surface Si is bonded to three O atoms with one O atom bonded to Cu. However, an alternative Cu/ C interface with interm ediate oxygen concentration can be constructed, which is done by removing two Cu bonding O atoms from half of the surface Si atoms and no t removing any oxygen atoms from the other half of Si surface atoms. The resulting Cu/ C interface is illustrated in Figure 4 3, which can be considered as a hybrid of Cu/ C:OO and Cu/ C:Si interfaces. The calculated adhesion energy of this hybrid Cu/ C interface from COMB potentials is 2.122 J/m 2 which is slightly smaller than the average adhesion energy of Cu/ C:O O and Cu/ C:Si interfaces, which is 2.227 J/m 2 This smaller than average adhesion energy results from larger Cu Si distances at the interface ; in other words this hybrid Cu/ C interface shows weaker Cu Si interactions than that of the Cu/ C:Si interfac e. Compared to the Cu/ C:O interface, on the other hand, this alternative intermediate
80 Cu/ C interface has a larger adhesion energy. This is indicative that this intermediate oxygen concentration hybrid interface is more energetically favorable than the original Cu/ C:O interface, and that oxygen atoms tend to saturate the bonding between Cu and Si atoms at the Cu/SiO 2 interface, rather than distribute themselves evenly throughout the interface. 4. 2 quartz Interfaces quartz is the most stable SiO 2 polymorph and amorphous silica is most widely used in electronic devices, we also investigate the properties of fcc Q (010) and fcc Cu (001)/a SiO 2 interfaces and predict the adhesion energies with the COMB potentials. We cr Q interfaces based on different interface terminations, OO O type I, O type II, and Si terminations. The OO terminated interface has surface O atoms bonded to Si and Cu atoms. Two types of O terminated interfaces are considered : type I and type II, which are constructed by removing one of the two Wyckoff positions occupied by O atoms. The Si terminated interface has half surface Si atoms bonded to the Cu slab, while the other half of the Si atoms are dangling. Force and energy m inimization are carried out to find the lowest energy interfacial configurations for the Cu/ Q:OO, Cu/ Q:O I, Cu/ Q:O II and Cu/ Q:Si interfaces, as shown in Figs. 4 4a and 4 4b. The minimized Cu O and Cu Si bond lengths at the interfaces from our MD simulations are ~1.91 and ~2.55 respectively, approximately 3 C interfaces. The predicted adhesion e nergies are given in Table 4 3. The longer Cu O and Cu Si bond lengths are reflective of the slightly smaller adhesion energies for the Cu/ Q:OO and Cu/ Q:Si interfaces compared to the Cu/ C:OO and Cu/ C:Si interfaces. However, both type I and type II Cu/ Q:O interfaces exhibit ~ 30% stronger adhesion than the
81 Cu/ C:O interfaces. The relatively large difference in adhesion energy for the Cu/ Q:O interfaces originated from larger surface energies for the O terminated Q surfaces, averaged 2.6 J/m 2 c ompared to that for the O terminated C surface, 1.3 J/m 2 Larger surface energies for the O terminated Q surfaces indicate less stable surfaces and also smaller (less negative) numbers for the term used in Equation 4 1 to calcul ate adhesion energies. This is the origin of the larger adhesion energies for the Cu/ Q:O interfaces. Overall, the trend for the adhesion energies is similar to the trend for the Cu/ C interfaces, and is consistent with the intuition that adhesion energy should increase with the number of oxygen atoms at the interface. 4. 2 .3 Adhesion of Cu/a silica Interfaces The bulk amorphous silica used in the Cu/a SiO 2 interfaces is prepared using the melt and quench method, as discussed in Shan [ 16 ] and Chapter 3. The a SiO 2 slab is created from the bulk a SiO 2 equilibrated at 300K for 50ps, and then combined with the Cu layer, followed by another 50ps of equilibration at 300K. The equilibrated structure of the Cu/a SiO 2 interface is illustrated in Figure 4 5. Examining the equilibrated interfacial structure, we note that the interface Cu atoms are oxidized and acquire charge values from 0.3 e 0.6 e. The oxidized Cu atoms and interface O atoms form a thin copper oxide layer with characteristics th at resemble the Cu 2 O phase. The atomic positions of the second Cu layer are only slightly influenced by the O atoms and the oxidized first Cu layer, and the charge values range from 0.0 e 0.1 e. From the third Cu layer inward, the Cu atoms retain their b ulk fcc positions and are charge neutral. This indicates that the Cu films oxidize in contact with glass substrates but that the oxidation is limited to the first two Cu layers. Examining the bonding between the Cu and a SiO 2 thin films, we
82 note that the a verage Cu O and Cu Si bond lengths are ~2.04 and ~2.15 respectively. We further predict that 22% of the interface Cu atoms are bonded to the glass substrate through Cu O bonds, and 16% are bonded through Cu Si bonds; the rest ( 62% ) of the interface Cu atoms are not bonded to the a SiO 2 These bonding statistics reflect the fact that O atoms may play more important roles than Si atoms in terms of bonding the Cu thin film to glass substrates. To investigate the effect of interfacial defects on the adhesi on energy of the Cu/a SiO 2 interface, we introduce oxygen vacancies to the interface and create three Cu/a SiO 2 interfaces with different numbers of interfacial oxygen defects (V O ) : 0, 10 and 20, respectively; 10 and 20 oxygen vacancies at the Cu/a SiO 2 in terface correspond to concentrations of 0.565 and 1.130 V O /nm 2 respectively. Force and energy are then minimized with the COMB potentials to calculate the adhesion energies, which are given in Table 4 4, and are compared to the values from experiments. 112, 113 The predicted adhesion energies of the Cu/a SiO 2 interface are reasonably comparable to the experimental values and show the same qualitative trends compared to Cu/ C and Cu/ Q interfaces Specifically, the adhesion e nergies decrease significantly with decreased concentration of oxygen atoms at the interface. As indicated in Table 4 4, the percentage of Cu O bonds at the Cu/a SiO 2 interface significantly decrease with the introduction of oxygen vacancies (with 20 vacan cies the Cu O percentage dropped to 11%), while that of the Cu Si bonds slightly increase. This bonding analysis indicates that Cu O bonds are the major contributor to the adhesion strength of Cu and SiO 2 substrates. This is due to the strong Coulombic att raction between oxidized Cu atoms and negatively charged O atoms. Additionally, although Cu Si bonds form at the
83 interface, the ionic repulsion strongly reduces their influence on adhesive strengths. Table 4 5 summarizes the average Cu O and Cu Si bond len gths at the Cu/SiO 2 interfaces for all three types of SiO 2 slabs considered 4. 3 quartz Interfaces To examine charge transfer across the interface and the effect of interface termination on the results we calculate the charge difference, along the interface normal direction Here, and are the distribution of charges of all of the atoms of the Cu/SiO 2 interface system, the individual Cu layer, and the individual SiO 2 slab, respectively. The specific Cu/SiO 2 interfaces investigated are the Cu/ C and Cu/ Q interfaces. The SiO 2 slab is removed from a fully relaxed, charge equilibrated Cu/SiO 2 interface and a charge equilibration with Cu atoms fixed is carried out to obtain The opposite process is carried out to obtain In other words, the physical meaning of the quantity indicates the ch ange in the charge quartz slabs an analogous analysis is found in Nagao [ 108 ] The charge differences, C interfaces are illustrated in Fig ure 4 6, while the Q:Si interfaces are shown in Figure 4 Q interfaces exhibit similar behaviors: the difference at the OO terminated interface is larger t han that at the other two interfaces, indicating the greatest charge transfer at the Cu/SiO 2 interfaces, particularly between Cu and O atoms. The interface with the largest number of Cu O bonds has the strongest adhesion C result is qualitatively consistent with electron density difference analysis found by DFT. 108
84 4. 4 Summary To conclude, we have applied empirical, variable charge COMB potentials for Si/SiO 2 and Cu/Cu 2 O quartz a nd Cu/amorphous silica interfaces with different types of terminations that mimic oxygen rich conditions intermediate and oxygen lean conditions, or different densities of interfacial oxygen defects. We also carried out DFT calculations, with both the LDA and the GGA, to cristobalite interfaces; results from these calculations were used to train the COMB potential mixing parameters. The results indicate that the COMB potentials reproduce the structural, adhesive, a nd electronic cristobalite interfaces reasonably consistently with DFT quartz interfaces are consistent with the notion that adhesion is strongest for oxygen rich conditions. The COMB potential was further applied to predict the adhesion energies of Cu/a SiO 2 interfaces with different concentrations of interfacial oxygen vacancies. The results indicate that the adhesion energy decreases as the number of oxygen vaca ncies increases due to decreased Cu O bonds The COMB potential for Si/SiO 2 /Cu/Cu 2 O is shown to be adequate for carrying out variable charge MD simulations of Cu/SiO 2 interfaces, and this should pave the way for significant progress in modeling these types of challenging interfaces.
85 Table 4 1. Properties of Cu 2 O calculated with the developed COMB potentials compared to that from experiments and DFT calculations. Properties Expt. 114 DFT PBE COMB a 0 () 4.274 4.267 4.226 E c (eV/Cu 2 O) 11.34 11.91 C 11 (GPa) 123 126 105 C 12 (GPa) 107 106 89 C 44 (GPa) 12 15 71 B (Gpa) 112 113 94 G (Gpa) 8 10 45 q Cu (e) 0.55 0.79 Table 4 2. Adhesion energies (in units of J/m 2 ) of the fcc cristobalite (001) interfaces from Ref. [ K. Nagao, 108 ] (first DFT LDA column), our DFT calculations (second DFT LDA column and DFT GGA column), and MD simulations using COMB potentials (COMB column). Cu/ C:OO, Cu/ C:O and Cu/ C:Si denotes OO O and Si term inated Cu/ C interfaces, respectively. Type of interface DFT LDA 108 DFT LDA DFT GGA COMB Cu/ C:Si 1.406 1.432 1.034 0.864 Cu/ C:O 1.555 1.591 1.222 1.734 Cu/ C:OO 3.805 3.987 3.601 3.591 Table 4 3. Adhesion energies (in units of J/m 2 ) of the fcc quartz (010) interfaces from MD simulations using COMB potentials. Cu/ Q:Si, Cu/ Q:O I, Cu/ Q:O II and Cu/ Q:OO denotes Si O type I, O type II and OO terminated Cu/ Q interfaces, respectively. Type of interface W (J/m 2 ) Cu/ Q:Si 0.850 Cu/ Q:O I 2.494 Cu/ Q:O II 2.646 Cu/ Q:OO 3.472 Table 4 4 Adhesion energies (in units of J/m 2 ) of the Cu/a SiO 2 interfaces with different numbers of oxygen vacancy defects from MD simulations using COMB potentials. Experimental adhesion energy values from the literature for Cu/a SiO 2 interface and the bonding analysis between Cu and a SiO 2 are also given. Type of i nterface W (J/m 2 ) Interfacial Cu bonding (%) Exp COMB Cu O Cu Si Non bonded Cu Cu/a SiO 2 + 0 V O 0.5 1.2 113 0.6 1.4 112 1.810 22 16 62 Cu/a SiO 2 + 10 V O 0.629 13 18 79 Cu/a SiO 2 + 20 V O 0.289 11 19 80
86 Table 4 5 Average Cu O and Cu Si bond lengths at the Cu/SiO 2 interfaces after equilibrated with COMB potentials. Type of interface Cu O bond lengths () Cu Si bond lengths () Cu/ C 1.86 2.32 Cu/ Q 1.91 2.55 Cu/a SiO 2 2.04 2.15
87 A B Figure 4 1. Relaxed structures of the fcc cristobalite (001) interfaces, OO O and Si terminated, respectively, with DFT PBE calculations viewing along A )  and B )  directions. Cu layer on the top, SiO 2 layer on the bottom (larger spheres are Si, smaller ones are O).
88 A B Figure 4 2. Minimized structures of the fcc cristobalite (001) interfaces, OO O and Si terminated, respectively, with MD simulati ons using the COMB potentials viewing along A )  and B )  directions. Cu layer on the top; SiO 2 layer on the bottom. Atoms are color coded by the charge values indicated by the color bar (navy is negatively charged, cyan is charge neutral, and red is positively charged).
89 A B Figure 4 3. Minimized structures of an alternative intermediate oxygen concentration fcc cristobalite (001) interface with MD simulations using the COMB potentials viewing along A )  and B )  directions. Color code is the same as that of Figure 4 2.
90 A B Figure 4 4 Minimized structures of the fcc quartz (010) interfaces, OO O type I, O type II, and Si terminated, respectively, with MD simulations using the COMB potentials viewing along A )  and B )  directions. Color code is the same as that of Figure 4 2.
91 Figure 4 5 M inimized structures of the fcc Cu (001)/a SiO 2 + 0 V O interfaces with MD simulations using the COMB potentials. Color code is the same as Figure 4 2.
92 Figure 4 6 Charge difference across the fcc cristobalite (001) interfaces. The Cu slab is to the right of the interface; the SiO 2 slab is to the left of the interface.
93 Figure 4 7 Charge difference across the fcc quartz (010) interfaces. The Cu slab is to the right of the interface; the SiO 2 slab is to the left of the interface.
94 CHAPTER 5 DEVELOPMENT OF COMB HF/HFO 2 POTENTIAL AND ITS AP PLICATION S TO AMORPHOUS HF O 2 AND SI/HFO 2 INTERFACES The continuous downsizing in the feature dimensions of microelectronic devices and ultra large integrated circuits (ULSIs) dimensions has necessitated ever thinner dielectric oxide layers. Silica ha s been used as the gate oxide for decades. However it had to become so thin (~1 nm) that the generation of leakage current had become a serious problem. To overcome this limitation, oxide materials with higher dielectric constants have replaced silica. Of all the oxide materials explored to date, hafnium dioxide (HfO 2 hafnia) is considered to be one of the most important because of its relatively high dielectric constants (typically 20 25 depending on the polymorphs), 115, 116 and high thermal stability when in contact with silicon (both semiconductor and poly Si gate el ectrode). 117 Moreover, it ha s been reported both experimentally and computationally that the incorporation of nitrogen 118, 119 and fluorine 120, 121 into hafnia films used in current devices further reduces the leakage current by passivating the oxygen vacancies. 122 However challenges still remain. In particular, the formation of a low dielectric constant hafnium silicate interfacial layer 123 and the presence of ox ygen defects at both the hafnia/silicon interface and in the hafnia film that trap electrons lead to a decrease of channel mobility and an instability in the threshold voltage. 124, 12 5 Hafnia crystallizes into three crystalline polymorphs at ambient pressure: a monoclinic phase, a tetragonal phase, and a cubic fluorite phase. The monoclinic to tetragonal transition takes place at 2000 K, while the tetragonal to cubic transition occurs at 2900 K. The melting point of the cubic phase is at 3085 K. 126 There are also two high pressure phases: the orthorhombic I phase ( Pbca,
95 Brookite type struc ture) above 10 GPa and the orthorhombic II phase ( Pnma PbCl 2 type or cotunnite structure) above 30 GPa. 127 Monoclinic hafnia has a band gap of 5.68 eV. 128 When it comes to empirical potentials for MD simulations, currently no published p otential predicts monoclinic hafnia to be the most stable phase or correctly describes the properties of hafnia polymorphs. In this work therefore, we develop a COMB potential that, by design, gives the correct phase order and reasonably well reproduces key properties of b oth elemental hafnium and hafnia polymorphs. Moreover, the potential is compatible with the previously developed potentials for the Si/SiO 2 system and amorphous silica, 16 allowing the simulation of HfO 2 /Si interfaces. This work would have not been made possible were it not for the daily discussions with Bryce Devine. Travis Kemper assisted in first principles DFT calculations of the elastic constants of the tetragonal hafnia phase. A senior undergraduate student, A. Dan ielle Holton ran MD simulations with the developed Hf/HfO 2 potential to obtain amorphous hafnia (a HfO 2 ) and characterized its properties. A first year graduate student, Xuan Sun, constructed various Si/a HfO 2 systems by combining the a HfO 2 structure with Si slabs and is currently conducting nanoindentation simulations of the Si/a HfO 2 systems. 5. 1 Potential Development and Properties of Hf Metal and H fO 2 Phases 5. 1 .1 Formalism and Fitting Procedures To model Hf and HfO 2 with the COMB potential formalism describe d in full details in Section 2.2 several correction terms were added to the main COMB formation which are described in Sections 2.2. 4 with Equations 2 4 2, 2 4 4, 2 4 5 2 4 6 and 2 4 7. We applied the multi phase fitting method, discussed in Section 2 2, to determine the
96 parameters. A total of 8 phases were included in the fitting procedure for the COMB potential for Hf/HfO 2 : the ambient pressure monoclinic, tetragonal, and cubic phases, hig h pressure orthorhombic I (OI) and orthorhombic II (OII) phases, rutile (P4 2 /mnm) PbO 2 (Pbcn) phases, and the pure Hf hexagonal metal phase (P6 3 /mmc). The Hf metal in fcc phase was not included in the fit because the fcc hcp energy difference is sole ly dependent upon the Legendre polynomial and can be incorporated after all the PbO 2 type phases do not appear in the phase diagram of hafnia, but some empirical potentials for AX 2 system have predicted t hem to be low energy states, or even the ground state. This has been an endemic problem in developing potentials for monoclinic zirconia, ZrO 2 129 131 Hence, we include these two phases in the fit to ensure that they are much higher in energy than the experimentally observed more stable phases. Since the parameterization of the Iwasaki potential for Si/HfO 2 was determined by fitting the atomic forces obtained with the potential to thos e obtained from first principles calculations on Ge/ZrO 2 or Si/HfO 2 interfacial models 28, 29 the potentials did not necessarily reproduce the structural and mechanical properties of the hafnium metal and hafnium oxide polymorphs as discussed in more detail below. However, the parameter set from the Iwasaki potential is still a va luable guideline and is used as the initial guess in our fitting procedure. The detemined parameterization for the Hf/HfO 2 COMB potential is given in Table A 1 The mixing rule scaling factors for Hf X where X can be any element, are set to 1.0. 5. 1 .2 Properties of H f Metal Table 5 1 compares the properties of hafnium metal predicted by the COMB potential with values from experiments, first principles calculations, and the Iwasaki
97 potential. Since Iwasaki did n ot publish these values in Iwa saki [ 28, 29 ] we implemented his potential in our in hous e HELL code and calculated the properties with the published parameters. The COMB potential not only predicts the lattice constants of Hf better than the Iwasaki potential, but it also yields almost exact agreement with experimental values for the elastic constants and the bulk and shear moduli. Most importantly, the COMB potential gives the correct order among the different phases. In particular for metallic Hf, the energy difference between the hcp and fcc phases is very close to that calculated from DFT calculations. This is actually straightforward to achieve with the COMB potential since fcc and hcp phases cannot be distinguished with first nearest neighbor potentials without additional correction terms. Fitting the coefficient for the Legendre polynomial of the COMB potential gives the right hcp to fcc phase order with the exact energy difference (as discussed in Section 2.2. 4.2) The Iwasaki potential does not contain a correction term such as that used here, and he nce predicts zero energy difference between the fcc and hcp phases. With an almost exact fit to the hcp to fcc energy difference, the COMB potential also predicts the unstable stacking fault energy as 1.695 J/m 2 which is very close to the value of 1.662 J/m 2 calculated from DFT. Also shown in Table 5 1 are linear thermal expansion coefficients, 6.910 6 K 1 and 5.910 6 K 1 predicted from the COMB potential and from experiments, respectively, and the defect formation energies of Hf vacancy and interstiti al. The COMB potential predicts the right trend and an interstitial formation energy that is very close to the DFT value; however, the COMB potential predicts a vacancy formation energy in Hf metal that is larger than that from DFT calculation by almost a factor of 2. This is also the case for the surface energies COMB predicts the
98 right trend for the surface energies but the values are approximately twice as large as the values from the DFT calculations for some of the surfaces considered. These overesti mations are consistent with the overestimation of the bulk cohesive energy. 5. 1 .3 Properties of H fO 2 Phases Table 5 2 compares of the properties of the monoclinic, tetragonal and cubic hafnia phases predicted by the COMB potential with experiment, DFT calc ulations and the Iwasaki potential. The properties from the Iwasaki potential were again calculated with our implementation of the potential. Using a damped force minimization (quench) method with the Iwasaki potential, the monoclinic phase spontaneously t ransforms to a PbO 2 type structure with a much lower cohesive energy than the monoclinic phase. As a result, the properties for the monoclinic phase predicted from the Iwasaki potential are left blank in Table 5 2 The Iwasaki potential predicts higher cohesive energies for the monoclinic and tetragonal phases than for the PbO 2 type structure; thus these phases are metastable in the Iwasaki potential. Not shown in the table are the relaxed surface energies for the (111), (001) and (100) surfaces of the monoclinic phase: 1.07, 1.17 and 2.02 J/m 2 respectively, compared to 1.25, 1.45, and 1.79 J/m 2 from a first principles study. 132 By construction, the COMB potential reproduces the phase order and elastic properties of the hafnia phases reasonably well. The monoclinic phase is the most stable structure with the COMB potential, which is in agreement with both experim ents and DFT calculations. This is demonstrated in Figure 5 1 (a), which illustrates the cohesive energies as a function of unit volume (E V curves) for various hafnia phases calculated from DFT with the Perdew Burke Ernzerhof exchange correlation function al 54 (DFT PBE), and Figure 5 1 (b), which shows the corresponding data for the
99 COMB potential. As is the case in the DFT calculations, the high pressure OI phas e has a lower energy than the tetragonal phase, and the equilibrium volume is substantially smaller than that of the monoclinic and tetragonal phases. The transition barrier should be high enough to prevent the undesired monoclinic to OI phase transition. Table 5 3 quantifies this comparison by comparing the cohesive energies of monoclinic, tetragonal, cubic, OI and OII hafnia phases calculated from DFT PBE and the COMB potential relative to the monoclinic phase. The DFT calculations of these known phases f rom this work are in excellent agreement with the literature, 133 and the correct phase order of the phases is again captured by the COMB potential. Also given in the table are PbO 2 PbO 2 type HfO 2 structure is predicted by the DFT calculations to have the second lowest PbO 2 type phase in ex periment may be attributed to the fact that large negative pressure is required for a monoclinic phase to transform to the PbO 2 type, as indicated by the PbO 2 system illustrated in Figure 5 1 (a). The fact that this pha se is predicted to be so low in energy could also be a short coming of the DFT PbO 2 type HfO 2 structures predicted from the COMB potential is opposite to that from the DFT PBE calculations, in both c ases their energies are sufficiently higher than the monoclinic phase to prevent their occurrences in MD simulations at standard temperatures and pressures For a PbO 2 type structures, negative pressure must be app lied (which corresponds to a volume expansion), accompanied by a decrease of the coordination number from 7 to 6 on all Hf atoms and 3/4 to 3 on all O atoms. The
100 monoclinic to PbO 2 transition was never observed in considered test cases, as discuss ed later in this section. The finite temperature properties are also of significant interest. Table 5 4 gives the thermal expansion coefficients of the monoclinic phase hafnia calculated with the COMB potential, compared with those from experiments. 134 The coefficients of thermal expansion in the a and c axes are slightly overestimated relative to the experimental values, while that in the b axis is in fairly good agreement. As a result, the predicted volumetric thermal expansion coefficient is ~30% larger than experimental values. We also tested the temperature dependence of the lattice parameters and axial angles of the monoclinic hafnia on heating and cooling using the COMB potential for Hf/HfO 2 developed in this work. The monoclinic phase was heated from 300 K up to 3500 K within 0.1 ns with three dimensional periodic boundary conditions applied the lattice parameters gradually increase due to thermal expansion and the axial angles remain fairly constant. No phase transformation was predicted while the monoclinic phase melted at around 3600 K. The lattice parameters and axial angles obtained from cooling the structure that was equilibrat ed at 3500 K follow the heating curves back down to 300 K. Since this system does not have any nucleation site for heterogeneous nucleation, this is not compelling evidence for the absence of a structural phase transition. To provide such a heterogeneous n ucleation site, surfaces were created on the monoclinic phase by lifting the periodic boundary condition on the (100) direction. This system, in which the surfaces act a nucleation sites, transforms to the tetragonal phase upon heating after 2.2 ps at 2300 K. On undergoing the phase transition, the 7 fold Hf
1 01 atoms become fully 8 fold coordinated, while the 3 fold O atoms become 4 fold coordinated. This is characteristic of a transition from the monoclinic phase to a fluorite like phase. This transition is q uantified by the change in time averaged lattice parameters in a, b and c axes upon heating from 2000 K and 300 0 K as illustrated in Figure 5 2. A t 2300 K (around 8 ps) the a axis spacing decreased from 5.236 to 5.196 while the b axis also decreased fr om 5.242 to 5.181 which corresponds to the lattice parameter of the tetragonal phase in the a b axis, indicating the phase transition (the measuring of a axis spacing takes the center most two unit cells to exclude the surface effect). This phase tra nsition is confirmed by coordination number change and the radial distribution function. No further phase transition was found upon further heating. A similar surface/interface enhanced tetragonal to monoclinic phase transition phenomenon has been shown in a first principles study by Carter et al. Al 2 O 3 /t ZrO 2 (001) interface. 135 5.1. 4 Defect Formation Energies in Monoclinic H fO 2 Table 5 5 compares the formation energies for the cation and anion vacancies, interstitials and Frenkel pairs and the Schottky defects in the monoclinic phase of HfO 2 determined from the COMB potential with DFT results. 136 The defect formation energies were calculated with the equation where is and are the energies of defective and perfect bulk structures, respectively, is the number of defects and is the chemical potential of the defect. The chemical potential of a hafnium atom in both the metallic and oxi de phases at 300 K was used in the equation to estimate hafnium containing defects, and that of O 2 molecule at the same temperature was used for defects that contain oxygen. The charge neutrality of the
102 system is maintained after creating the defects. For example, if we consider a Hf vacancy, a Hf ion, instead of a Hf atom, is deleted from the structure and it leaves its charges (+3.48 e) in the system. In contrast, when we create a Hf interstitial, a neutral Hf atom is added to the system. A cation Frenkel defect is created by introducing both a n Hf ion vacancy and a Hf ion interstitial. Due to the dynamic charge transfer feature inherent to the COMB potential, created defects are able to take on charges from the surrounding ions while maintaining charge n eutrality over the system as a whole. Charges associated with interstitial defects after charge equilibration are given in parenthesis in Table 5 5 Two values are given for oxygen interstitials, indicating the incorporation of oxygen interstitials into 3 or 4 fold coordinated oxygen sites, respectively. Obtaining a good estimate of defect formation energies of the oxygen defects is particularly important since they act as charge traps and diffuse in the oxide as interstitials. Since the defect formation energies were not part of the fitting set of physical properties, the fact that the COMB potential predicts the defect formation energies with fidelity similar to that of DFT can be viewed as an indication that much of the correct physics is being captured Removing a hafnium atom from the bulk monoclinic phase results in the breaking of 7 Hf O bonds; thus the energy required to form this point defect is the greatest among Hf/O vacancies and interstitials. Similarly, creating an oxygen vacancy requires brea king 3 to 4 Hf O bonds; hence the formation energy is less than that of the hafnium vacancy. Hf and O interstitials require less energy to form since increasing the coordination number is more energetically favorable than breaking several Hf O bonds. Of al l the point defects, the oxygen interstitial has the smallest formation energy. Also
103 shown are the estimated defect formation energies of Schottky defects, and cation and anion Frenkel pairs. T he energy required to create a Schottky defect is slightly larg er than that of one Hf and two O vacancies added together. Likewise, the formation energy of the cation Frenkel pair is little larger than the sum of one Hf vacancy and one Hf interstitial defect, which indicates that the association energy of the Hf vacan cy and interstitial defects are positive; this implies that these defects like to be far apart and additional energy is required to bring them close to each other. In contrast, the formation energy of the anion Frenkel is slightly smaller than the sum of c onstituent point defects, indicating slightly negative association energy. 5. 1 .5 Properties of Amorphous H fO 2 There are various techniques for depositing hafnia films including atomic layer deposition (ALD), 137 pulse laser deposition (PLD), 138 and chemical vapor deposition (CVD). 139 Deposited hafnia thin films are typically amorphous as deposited and polycrystalline after a post depo sition annealing. 132 It has been reported that hafnia films deposited at 300 C by ALD are predominantly monoclinic with texture. 140, 141 Due to their anisotropic characteristics, amorphous hafnia thin films deposited on Si substrates are more desired than the crystalline forms. It is therefore important for the developed COMB Hf/HfO 2 potential to be able to obtain and describe amorphous hafnia. Working with Danielle, w e apply the COMB Hf/HfO 2 potential to obtain amorphous hafnia via the melt and quench method similar to that described in Section 3.3.1. The hafnia hot melt obtained at 5000K was quenched to 0K and equilibrated at 300K. The radial distribution function of the obtained amorphous hafnia at 30 0K is shown in Figure
104 5 3. The major peak at 2.11 corresponds to primary Hf O bonds, and the peak at corresponds to the secondary Hf O bonds. The minor peaks at 3.19 and 3.71 correspond to O O and Hf Hf bonds, respectively. More work is currently ongoing in collaboration with Danielle Holton, to investigate the dielectric properties and diffusion coefficients of amorphous hafnia at various temperatures. 5. 2 Applications to Si/HfO 2 Interfaces Although hafnia is an excellent material for the gate ox ide dielectric, problems still exist; in particular the diffusion of oxygen through the hafnia layer to form a silica layer at the HfO 2 /Si interface. 142, 143 The formation of this interfacial layer between hafnia and silicon substrates, which decreases th e dielectric constant and leads to material deterioration, must be suppressed. In this section, we demonstrate the ability of the COMB potential to model the HfO 2 /Si interface and oxygen transport at the interface. To carry out this study, the Hf/HfO 2 pote ntial is implemented together with the COMB potential for Si/SiO 2 16 This is possible because the functional form and the oxygen parameters are exactly the same between the two potentials, thus allowing the oxygen to interact with both elements in a well def ined manner within the system. The short range interaction between Hf and Si follows the Lorentz Berthelot mixing rule described in Equations 2 2 1 to 2 2 4 and the Hf and Si ions also interact via long range Coulombic interactions. For this proof of principle demonstration of the potential developed in this work, the oxygen terminated cubic HfO 2 (100)/Si (100) interface and moclinic HfO 2 (100)/Si(100) interface were chosen.
105 5. 2 .1 Si and Cubic HfO 2 Interface 5. 2 .1. 1 Structural propert ies of the interface Figure 5 4 illustrates snapshots of the c HfO 2 :O/Si interface at (a) its initial condition (since the lattice parameters of c HfO 2 and Si are fairly close, the interface is fully coherent with 8.1% strain applied to the Si substrate i n the x and y directions), and (b) after 2 ps of evolution at 300 K. Initially the c HfO 2 layer carried +3.30 e on the Hf atoms and 1.485 e on the O atoms. The Si substrate is neutral and the system as a whole is charge neutral. T he interfacial Si and O atoms from the Si substrate and c HfO 2 :O, respectively, move slightly into the interface, forming Si O bonds. As a result a very thin SiO 2 layer with oxidized Si atoms is formed. The driving mechanism for the formation of the SiO 2 layer may be a result of the unfavorable coordination of Si atoms at the interface, which are 6 fold coordinated and have relatively high energy. Once the small displacement between the interfacial Si and O atoms takes place to form the SiO 2 layer, the interfacial Si atoms are 4 f old coordinated, the second layer Si atoms are 3 fold coordinated, the interfacial Hf atoms are 6 fold coordinated, and the oxygen atoms are 3 or 4 fold coordinated. This results in a more stable bonding environment at the interface than the starting stru cture. Figure 5 5 illustrates the distribution of charge at the c HfO 2 /Si interface for Si, O, and Hf atoms after 2 ps of evolution. At the interface, the change in charge on Si atoms is +2.308 e on average, while that for O and Hf atoms are 0.413 e and +0.383 e, respectively. The results indicate that interfacial Si atoms oxidize, forming Si O bonds and donating electrons to the O atoms in the cubic hafnia. The oxidation of Si atoms at the interface to form a silica layer is consi stent with experimental results; 142, 143 as we shall see, these strong Si O bonds contribute to the large value of the work of adhesion.
106 5. 2 .1. 2 Adhesion of the interface The work of adhesion quantifies the energy saving from two surfaces forming an interface, and is an indication of the stability and strength of the interface. The larger the work of adhesion of an interface, the more energy is required to separate the inte rface into two surfaces. Here the work of adhesion W of the c HfO 2 /Si interface was calculated from the standard equation : where and are the total potential energy of the relaxed fixed vol ume HfO 2 and Si slabs, respectively, is the energy of the relaxed interfacial structure, and A is the surface/interface area. The same method was applied to calculate the work of adhesion of Cu/SiO 2 interfaces with DFT. 108 All structures are equilibrated at 300 K, prior to quenching at T=0 K. The calculated work of adhesion of the quenched c HfO 2 /Si is 7.36 J/m 2 which indicates a relatively strong interface, and is in excellent agreement with that calculated from DFT PBE, 7.19 J/m 2 The surface energies calculated from the relaxed c HfO 2 (100) and Si (100) slabs with COMB potential are 5.74 and 3.03 J/m 2 respectively, compared to those obtained from DFT calculations, which are 5.56 and 2.83 J/m 2 respectively. Note that the Si (100) surface energy is high because of the 8.1% strain applied to (010) and (001) directions. The agreement between the COMB prediction and the DFT calculation indicates that the combination of COMB po tentials for Hf/HfO 2 and Si/SiO 2 is adequate for modeling Si/HfO 2 interfaces. However, additional extensive testing is required to identify cases for which this approach works less well.
107 5. 2 2 Si and Monoclinic H fO 2 Interface Since monoclinic phase is t he most stable hafnia polymorph and post deposition annealing of HfO 2 deposited on Si are predominantly the monoclinic (m HfO 2 ) phase 141, 144 we apply the COMB potentials to examine the properties of Si (100) /m HfO 2 (100) interfaces. A snapshot of such a system equilibrated with COMB po tentials is given in Figure 5 6 Interfacial Si atoms are oxidized, forming a thin SiO 2 layer which is consistent with experimental findings This preliminary result indicates the ability of COMB potentials to investigate a wide range of Si and HfO 2 (cubi c, monoclinic or amorphous) interfaces and systems. 5. 2 3 Nanoindentation of Si and Amorphous HfO 2 Systems HfO 2 is being applied in microelectronic devices; therefore in addition to electrical properties of Si/HfO 2 interfaces, mechanical properties of such interfaces are also important. As a result, we apply COMB potentials to investigate the mechanical properties of Si/a HfO 2 interface via nanoindentation simulations in a similar manner described in Section 3 6. Wo rking with Xuan, we took previously determined amorphous hafnia structure and made two a HfO 2 slabs with different thicknesses, 12 and 18 respectively. The equilibrated a HfO 2 slabs are combined with 27 thick Si slabs, followed by constant volume eq uilibration at 300K. Equilibrated Si/a HfO 2 interfacial structures are shown in Figure 5 7 The structures are indented with rigid, hemi spherical Si indenters with 33 diameter of curvature at 300K and as a preliminary test, the loading and unloading ra te s are 20 m/s. The load displacement curves of the two tested Si/a HfO 2 structures are illustrated in Figure 5 8. The loading depths are approximately two thirds of the a HfO 2 thicknesses, so that different displacements are provided in the figure. The Si/a HfO 2 struct ures exhibited permanent damage from
108 plastic deformation, which are identified from the load displacement responses. The extracted hardness and reduce modulus of the Si/a HfO 2 structures are given in Table 5 6. The 12 thick a HfO 2 interface shows smaller values for hardness and reduced modulus, due to a thinner a HfO 2 slab. This work is currently being redone with slower loading/unloading rates at 1 m/s, as well as being replicated on Si/a HfO 2 structures with larger surface areas of Si/a HfO 2 and a HfO 2 /vacuum interfaces than the ones considered here 5. 3 Summary The COMB potential for Hf/HfO 2 developed in this work captures most of the physical properties of Hf metal and hafnia polymorphs. However, it is necessary to sacrifice accuracy in some properties, such as the Hf metal lattice constants and cohesive energy, in order to guarantee the correct phase order among hafnia polymorphs. Most importantly, however, t he COMB potentials for Hf/HfO 2 Cu/Cu 2 O 30 and Si/SiO 2 16 can be seamlessly coupled together in MD simulation studies of real device size systems such as HfO 2 /Si interfaces, Hf x Si 1 x O y films, hafnia films on Si or SiO 2 and the entire gate stack (Si/SiO 2 /HfO 2 /Cu 2 O/Cu) found in MOS/CMOS devices The potential is also capable of describing amorphous hafnia and Si/a HfO 2 interfaces as well as carrying out nanoindentation studies on such systems. Thus, the COMB potential for Hf/HfO 2 should prove to b e a useful tool and an effective method for carrying out large scale MD simulations and computational studies.
109 Table 5 1 Properties of hafnium metal given by the COMB potential for Hf/HfO 2 developed in this work in comparison with experimental, densi ty functional theory (DFT) calculations and the Iwasaki potential Properties Exp COMB DFT (PBE) Iwasaki a 0 ( ) 3.1950 3.167 3.1935 3.1389 c 0 ( ) 5.0542 5.147 5.0533 5.1276 E c (eV/atom) 6.99 6.98 9.96 6.59 C 11 (GPa) 190 195 201 C 12 (GPa) 75 75 69 C 13 (GPa) 66 65 63 C 33 (GPa) 204 209 220 C 44 (GPa) 60 53 58 C 66 (GPa) 58 61 66 B (GPa) 110 112 114.59 112 G (GPa) 61 60 65.1 E ( fcc hcp ) (eV/atom) 0.055 0.060 0.0 E ( bcc hcp ) (eV/atom) 0.188 0.155 2.273 E ( dia hcp ) (eV/atom) 3.98 2.07 3.768 Unstable stacking fault (J/m 2 ) 1.695 1.662 2 ) 2250 1133 3006 (mJ/m 2 ) 2466 2003 4002 (mJ/m 2 ) 2936 3044 3752 Vacancy (eV/defect) 4.62 2.55 Interstitial (eV/defect) 6.17 6.23 6 K 1 ) 5.9 6.9 Denotes fitted properties
110 Table 5 2 Fitted properties of monoclinic, tetragonal and cubic hafnia phases given by the COMB potential for Hf/HfO 2 developed in this work in comparison with that from experiments, first principles calculations, a Buckingham type potential and the Iwasaki potential Monoclinic HfO 2 Properties Expt. COMB DFT Iwasaki a 0 ( ) 5.12 a 5.13 5.12 b, c -b 0 ( ) 5.17 a 5.21 5.20 b, c -c 0 ( ) 5.30 a 5.11 5.28 b, c -(deg) 99.2 a 98.8 99.7 b, c -E c (eV/HfO 2 ) 30.89 30.56 -B (GPa) 235 251 d -G (GPa) 120 -q Hf (e) 3.48 3.60 c 3.38 Tetragonal HfO 2 Properties Expt. COMB DFT Iwasaki a 0 ( ) 5.151 a 5.03 5.059 c 4.921 c 0 ( ) 5.29 a 5.05 5.1996 c 6.26 E c (eV/HfO 2 ) 30.79 30.40 41.85 B (GPa) 210 e 281 f 307 197 c 240 585 G (GPa) 154 128 342 C 11 (GPa) 567 495 1151 C 12 (GPa) 179 152 302 C 13 (GPa) 176 115 302 C 33 (GPa) 564 397 1149 C 44 (GPa) 126 90 288 C 66 (GPa) 128 125 288 q Hf (e) 3.46 3.33 c 3.30 Cubic HfO 2 Properties Expt. COMB DFT Iwasaki a 0 ( ) 5.08 a 5.02 5.05 b 6.25 E c (eV/HfO 2 ) 30.76 30.32 40.74 B (GPa) 295 201 c 280 d 556 G (GPa) 200 229 g 324 C 11 (GPa) 561 578 g 1098 C 12 (GPa) 161 121 g 285 C 44 (GPa) 111 83 g 270 q Hf (e) 3.44 3.49 c 3.22 Properties for the monoclinic phase are left blank, see text for explanation a Ref 134, 145, 146 b Ref 147 c Ref 133 d Ref 148 e Ref 127 f Ref 149 g Ref 150
111 Table 5 3 Cohesive energy of monoclinic hafnia (in units of eV/HfO 2 unit) and that relative to monoclinic hafnia for six other polymorphs calculated from DFT PBE and COMB potential for Hf/HfO 2 developed in this work. (eV/HfO 2 ) Monoclinic PbO 2 OI Rutile Tetragonal Cubic OII DFT a 30.45 -+ 0.065 -+ 0.156 + 0.237 + 0.385 DFT PBE (This work) 30.56 + 0.058 + 0.070 + 0.109 + 0.159 + 0.238 + 0.386 COMB 30.69 + 0.076 + 0.085 + 0.067 + 0.178 + 0.293 + 0.430 a Ref 133 Table 5 4. Axial and volume thermal expansion coefficients of the monoclinic hafnia calculated with the COMB potential developed in this work com pared with those from experiments a 0 b 0 c 0 V Exp a 6.1 ~ 8.4 1.2 ~ 2.7 9.6 ~ 12.6 20.6 ~ 23.5 COMB 10.52 2.11 13.89 29.46 a Ref 134 Table 5 5. Defect formation energies (in units of eV) of the monoclinic HfO 2 phase predicted with the COMB potential developed in this work compared with those from DFT calculations. The chemical potential of Hf in both metal and oxide phases were used to estimate the formation energies of Hf vacancy, interstitial and Schottky defect. Values in parenthesis are the charges associated with interstitial defects from DFT and COMB potential after charge equilibration, respectively. Defect DFT a COMB metal oxide metal oxide Hf Vacancy 16.9 5.7 23.53 6.49 Hf Interstitial -9.78 (+1.983e) 26.81 (+1.983e) O Vacancy 9.34, 9.36 11.95, 13.44 O Interstitial 4.2, 5.8 (0.0e) 3.97, 5.24 ( 0.002e) Schottky -47.10 30.07 Cation Frenkel -43.29 (+1.991e) Anion Frenkel -18.06 (+0.101e) a Ref 136 Table 5 6. Hardness and reduced modulus of Si/a HfO 2 structures with 12 and 18 thick a HfO 2 layers obtained from nanoindentation simula tions using COMB potentials developed in this work. Courtesy of Xuan Sun. Type of interface H (GPa) E r (GPa) 1 2 Si/a HfO 2 33.3 333.6 18 Si/a HfO 2 35.8 458.0
112 A B Figure 5 1. Cohesive energies as a function of unit volume (E V curves) for various hafnia phases calculated from A) density functional t heory (DFT) at the level of PBE and B) that calculated from the COMB potential for Hf/HfO 2 developed in this work COMB
113 Figure 5 2. Temperature dependence of the lattice constants of monoclinic hafnia with (100) surface on heating from 2000 to 3000 K showing the monoclinic to tetragonal phase transition at 2300 K using the COMB potential for Hf/HfO 2 developed in this work
114 Figu re 5 3. Radial distribution function of amorphous hafnia obtained with the melt and quench method using the COMB Hf/HfO 2 potential developed in this work. Courtesy of Danielle Holton.
115 A B Figure 5 4 Snapshots of the cubic HfO 2 and Si (c HfO 2 /Si) interfacial structure evolving at 300K (a) 0 ps, (b) 2.0 ps. Large (red) atoms are Hf, medium sized atoms (cyan) are Si and small atoms (blue) are oxygen.
116 Figure 5 5 Charge distribution of the c HfO 2 /Si int erface for Si (square symbols), O (circle symbols) and Hf (triangle symbols) atoms after 2 ps evolution. Length ()
117 Figure 5 6 Snapshots of the moclinic HfO 2 and Si (m HfO 2 /Si) interfacial structure evolving at 300K. Large (red) atoms are Hf, medium sized atoms (cyan) are Si and small atoms (blue) are oxygen.
118 A B Figure 5 7 Snapshots of Si/a HfO 2 systems with A ) 12 and B ) 18 thick a HfO 2 layers from nanoindentation simulations using COMB potentials a HfO 2 laye r on the top, Si slab on the bottom. Atoms are color coded by the atomic charges. Courtesy of Xuan Sun.
119 A B Figure 5 8 Load displacement curves of Si/a HfO 2 systems with (a) 12 and (b) 18 thick a HfO 2 layers from nanoindentation simulations. Loading and unloading rates are 20 m/s; loading curves shown in red, unloading curves shown in blue and green. Courtesy of Xuan Sun.
120 CHAPTER 6 DEVELOPMENT OF COMB TI/TIO 2 POTENTIAL AND ITS AP PLICATION S TO TIO 2 SUR FACES AND TIO 2 SUPPORTED COPPER CLUSTERS AND ISLANDS Titania (TiO 2 ) has attracted much attention for its diverse use in various industries and wide range of desirable surface properties. It has been utilized in such applications as pigments, gas sensors, s olar cells, photo catalysts and heterogeneous catalyst. 151 Titania is well known to crystallize into three natural polymorphs: rutile (space group: P4 2 /mnm ), anatase ( I4 1 / amd ) and brookite ( Pbca ). In addition to these naturally occurring phases, several high pressure phases are also found experimentally or sugges ted by theoretical calculations: columbite (PbO 2 space group Pbcn ), cotunnite (PbCl 2 Pnma ), baddeleyite (ZrO 2 P2 1 /c ) and fluorite (CaF 2 Fm 3m ). 152 A number of empirical, fixed charge rigid ion potentials for carrying atomistic and molecu lar dynamics simulations have been proposed for TiO 2 153, 154 Of these potential models, it i s considered that a formal charged model by Collins et. al 153 and a partial charged model by Matsui et al. 155 describe the relative phase stabilities of titania polymorphs and low Miller index rutile surfaces in a reasonable manner. 152 Since Rapp and Goddard developed the charge equilibration (QEq) method, 26 several empirical potentials adopting this variable charge scheme have also been developed for TiO 2 : MS Q, 152 SMB Q 156 and SMTB Q. 157 MS Q combines a short range two body Morse Stretch potential with a long range Coulombic, SMB Q combines the Coulombic with a short range second mome nt Buckingham potential, while the SMTB Q model utilizes a second moment tight binding scheme for the short range interactions. These potential formalisms have their strengths, however it is difficult for them to model both Ti metal and TiO 2 within the sam e functional form, mostly due to the significant difference in bonding between the two materials: metallic versus ionic. Due to this limitation, it is also
121 difficult to model metal clusters, such as Cu Al, Au and Ag interacting with TiO 2 with these poten tials. Therefore we propose a COMB potential for both Ti metal and TiO 2 polymorphs, and, that when combined with the Cu/Cu 2 O potential is able to model systems such as TiO 2 supported Cu metal clusters. This work received great assistance from Rakesh Behera. In particular, he supplied several titania polymorphs for the phase stability calculations, as well as characterized several empirical potentials in terms of reproducing titania properties. Patrick Chiu assisted in carrying out DFT calculations of defect formation energies at a grain boundary of rutile titania. A first year graduate student, Jackelyn Martinez, constructed the reconstructed and step rutile (110) surfaces and investigated their surface energies and relaxed geometries. 6. 1 Potential De velopment and Properties of Ti Metal and TiO 2 Phases 6. 1 .1 Formalism and Fitting Procedure The cosine bond bending terms, described in Section 2.2. 4.1, are applied to Ti O Ti and O Ti O bond angle s, and the additional repulsion term, described in Section 2.2. 4.3, is applied to Ti O bonds within the standard COMB formalism We apply the multi phase fitting algorithm described in Section 2.6 to determine the potential parameters. We have included a total of 6 phases in the fitting procedure for the COMB pote ntial for Ti/TiO 2 : the pure Ti hexagonal metal phase (P6 3 /mmc), rutile, anatase brookite cristobalite The Ti metal in fcc phase was not included in the training database because the fcc hcp energy difference is solely dependent upon the Legendre polynomial and can be incorporated after all the other parameters have cristobalite type phases do not appear in the phase diagram of titania, but for some parameter set s these phases have been predicted to be
122 low energy states, or even the ground state. Hence, we include these two phases in the fitting procedure to ensure that they are higher in energy than the experimentally observed phases. The final parameterization for the COMB Ti/TiO 2 potential is gi ven in Table s A 1 and A 2. 6.1. 2 Properties of Ti Metal Table 6 1 compares the fitted properties of hcp Ti metal from experiments 158 160 with the fitted and predicted values from first principles calculations 158 a MEAM potential 158 and the COMB potential (developed in this work). The lattice parameters from the COMB potential show 0.5% and +1.5% deviation from experimental values for a 0 and c 0 respectively. These deviations result in a larger c/a value of 1.620 compared to exper imental and DFT values. However, the cohesive energy and elastic constants are reproduced within 2% deviation with respect to DFT values. The fcc hcp energy difference is exactly fitted, and the predicted bcc hcp energy difference, 67 meV/atom, agrees well with DFT calculations and indicates the bcc phase is less stable than the fcc phase, while the hcp phase is the ground state phase. Point defects (vacancy and interstitial) and planar defects (stacking faults and free surfaces) are also predicted with th e COMB potential. The defect formation energy for a Ti vacancy is reasonably predicted while the Ti interstitial (octahedral site) is overestimated by ~ 55%. Stacking fault energies of hcp metals are important factors not only to their mechanical propertie s, but also to deformation mechanisms such as slips systems and dislocations. There are three basic types of stacking faults in hcp metals. Intrinsic type 1 ( I 1 ) changes the normal stacking sequence of to ABAB|CBCB ; type 2 ( I 2 ) changes it to ABAB|C ACA ; while an extrinsic stacking fault has a stacking sequence of ABAB|C|ABAB The unstable stacking fault (USF) energy is the
123 minimum energy barrier between the normal hcp stacking sequence and the I 2 stacking fault. These stacking fault energies of Ti ar e well described with the COMB potential; both I 2 and USF energies are within 13% of the DFT results. The surface energies of the three low Miller index surfaces are also given in Table 6 1. The COMB potential predicts a different trend for these surface energies compared to DFT values; however the predicted surface stability trend is the same as that from the MEAM potential. Predicted surface energies from the COMB potential are ~35% smaller than that from DFT calculations. 6. 1 .3 Properties of TiO 2 Phases and Phase Stability The properties of rutile, anatase and brookite phases from the COMB potential, compared to the experimental values, first principles calculations, a tight binding (SMTB Q) 1 57 method, and three other classical interatomic potentials (SMB Q 156 MS Q 152 and Matsui 155 ) are presented in Table 6 2. Our first principles calc ulations are carried out with the Vienna ab initio Program (VASP) within the generalized gradient approximation (GGA) 65 using the Perdew Burke Ernzerhof 54 (PBE) exchange correlation functional. We use plane wave basis sets with a 580 eV energy cutoff, a 444 Monkhorst Pack k point mesh 67 and projector augmented wave (PAW) ps eudopotentials 63, 64 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 4.0 10 4 1 for energies and forces, respectively. The Matsui and MS Q values are generated with the General Utility Lattice Program (GULP) 161, 162 using parameters provided in Matsui [ 155 ] and Swamy [ 152 ], respectively.
124 For the rutile phase, the COMB potential reproduces the lattice parameter, cohesive energy and elastic moduli with a 20% deviation compared to experimental and DFT values. Latt ice parameters, cohesive energies and elastic moduli of anatase and brookite phases from the COMB potential are also well reproduced, however the individual elastic constants of anatase and brookite, in some cases, show more significant deviation. For anat ase C 12 shows a ~ 60% deviation and C 44 shows a ~ +90% deviation. For brookite C 13 shows a ~ 30% deviation and C 66 shows a ~ +95% deviation relative to the DFT values. Compared to values from MS Q and Matsui potentials, the COMB TiO 2 potential reproduces worse elastic constants for anatase C 12 and C 44 and brookite C 13 but better brookite C 66 In fact, it can be seen that the shear elastic constants for all three phases from empirical potentials show substantial deviations from experimental and DFT values This may indicate the insufficiency of current empirical potential functions, including MS Q, Matsui and COMB potentials, to correctly reproduce shear elastic constants for TiO 2 polymorphs. One insufficiency of the potential function may be the angular f unction, and a new model for the angular function may be required to correctly reproduce TiO 2 shear moduli. The phase order of titania polymorphs, including naturally occurring phases, high pressure phases and several hypothetical AB 2 phases, from the COMB potential compared to first principles DFT calculations and other empirical potentials are provided in Table 6 3. COMB and SMB Q potentials predict the rutile phase to be the most stable phase, followed by anatase and brookite phases. MS Q and Matsui pote ntials, while predicting rutile as the ground state phase, predict that brookite is lower in energy than anatase. The COMB potential does not predict exactly the same phase order for all
125 other AB 2 phases compared to first principles calculations, however these high pressure or hypothetical AB 2 phases are correctly predicted to be much less stable than the rutile, anatase and brookite phases. 6.1 .4 Defect Formation Energies of Bulk Rutile TiO 2 Point defects govern the electronic conductivity and electrical properties of TiO 2 Point defect formation energies are therefore calculated with the developed COMB potential and compared to experimental values 163 and first principles calculations. 164 The calculated defects include Ti and O vacancies, Ti and O interstitials, cation and anion Frenkel and Schottky defects. For COMB potential calculations, system charge neutrality is maintained after creating defects. The calculated defect formation energies are compared to literature values and are presented in Table 6 5 In particular, defect formation energies from the COMB potential are provided for oxidized (Ti in oxide & O 2 gas), standard (Ti metal & O 2 gas), and reduced (Ti metal & O in oxide) conditions. The experimental values for the defect formation en ergies are their upper limits. 163 In general the COMB potential ove restimates the defect formation energies by factors from 2 to 5 compared to experimental values and first principles calculations. This is due to strong ionic contribution resulted from larger fractional charges. Under the standard condition, COMB predicts defect formation energies that have the same trend as the values predicted from first principles calculations. The calculated trends for Schottky and Frenkel defects also agree with the trends from experiment and DFT calculations. In particular, formation of a Schottky defect requires the most energy in bulk rutile, while the formation of an anion Frenkel requires the least energy. The results from COMB potential calculations also indicated that more energy is required to create Ti vacancies than O vacanc ies. This is because a Ti vacancy is
126 created from breaking six strong Ti O bonds, while an O vacancy requires breaking only three; hence the Ti vacancy formation energy is larger than of the value for an O vacancy. Ti or O interstitials, on the other hand, are not significantly affected by strong Coulombic interactions since interstitials attract and repel dissimilar and similar atoms, respectively. The O interstitial forms an almost linear Ti O Ti bond angle with the surrounding Ti atoms, so creating this defect requires additional energy to overcome the b ond bending term, defined in Equation 2 3 9. On the other hand, the Ti interstitial forms an O Ti O bond close to the ideal bond angle hence it is more energetically favored than the O interstitial. 6. 1 .4 D efect Formation Energies at Rutile TiO 2 Grain Boundaries Grain boundaries in polycrystalline materials are well known to affect their electronic properties, 165 167 and defect comp (210) grain boundaries significantly affect the electronic properties of rutile TiO 2 168 rutile TiO 2 using COMB potential, and compare the results to that of first principles DFT calculations. 169 The grain boundary structures are built following the method described in Chiu [ 169 ] and the results are presented in Table 6 5 The relaxed structures are also illustrated in Figure 6 1. Compared to DFT values, the formation energies for cation Frenkel and Schottky defects from COMB potential are overestimated by several factors. These over estimations in defect formation energies also originate from strong ionic contribution due to larger fractional charges. Nevertheless the cation Frenke l is predicted to be more favorable than Schottky defect, which is consistent with experimental values.
127 For all three conditions considered, COMB potential predicts the mixed A structure to be the most stable Schottky defect, followed by clustered and mixe d B configurations of Schottky defects, consistent with DFT calculations. The significant difference in formation energies among these Schottky defects from COMB potential can be attributed to the number of under coordinated Ti and O atoms more under coo rdinated atoms larger the formation energies. An analysis of the coordination number of Ti and O atoms of the system is presented in Table 6 6 The mixed A configuration Schottky defect contains 23.8% 5 fold coordinated Ti atoms, and the clustered Schottky defect contains 4.8% 4 fold and 19.0% 5 fold coordinated Ti atoms. This difference in under coordinated atoms results a more preferred mixed A than clustered Schottky defect. Comparing the mixed B and distributed Schottky defects, the percentage of atoms that are normally coordinated decreased to 77.0% and 78.6%, respectively. The significant increase in the numbers of under coordinated Ti and O atoms for the mixed B and distributed configurations of Schottky defects is responsible for their higher defect formation energies. Examining the Frenkel defects, the cation Frenkel pair is predicted to be more preferred than the anion Frenkel pair from COMB potential, which agrees with DFT calculations. This also can be attributed to the overall number of under coo rdinated Ti and O atoms of these two Frenkel types of defect systems: the number of under coordinated Ti and O atoms in the cation Frenkel is less than that of the anion Frenkel. Thus, the defect formation energies at rutile grain boundaries from COMB pote ntial are consistent with available experimental values and first principles DFT calculations.
128 6. 2 Surface Energies of TiO 2 Phases Single crystalline TiO 2 surfaces are important in wide range of applications 151 An important aspect of assessing the quality of a TiO 2 empirical potential is its ability to describe surfaces. Therefore, we investigated the quality and capability of the developed COMB Ti/TiO 2 potential by calculating the surface energies of rutile, reconstructed rutile (110) and anatase. Since experimenta l values for surface energies are difficult to obtain, we compare COMB potential values to our or literature first principles DFT calculation and other empirical potential values 6. 2 .1 Low Miller Index Rutile Surfa ces Table 6 7 compares the surface energi es of (110), (100), (101) and (001) surfaces of rutile phase predicted from the COMB potential, first principles calculations 170 174 and empirical po tentials. 155, 156 MS Q values are calculated with GULP using given parameters. 155 For un relaxed surface energies the COMB potential predicts the same trend compared to DFT calculations, in particular to Ramamoorthy et al. 174 The relaxed surface energies generally follow the same trend, except that the (100) surface is predicted to be more stable than the (110) surface after relaxation. Interestingly, the MS Q potential predicts an incorrect o rder for (100) and (001) surface energies. The SMB Q potential, on the other hand, does a good job in predicting the surface energies for all the low index surfaces considered here. Examining the relaxed surfaces from the COMB potential, we found the under coordinated surface Ti atom being pulled into the bulk, and this is also predicted for the (001) surface. This over relaxation of surface Ti atoms for (100) and (001) surfaces are responsible for the 60% and 53% decreases, respectively, between the un relaxed and relaxed surface energies. On the other hand, the decrease in surface energies for (110) and (101) surfaces are much smaller at 28%
129 and 42%, respectively. Although the relaxed (100) surface is predicted to be the most stable from COMB potential, the (110) surface is structurally stable during relaxation at higher temperatures and does not evolve to other surface configurations. Therefore the COMB potential should be adequate for modeling the rutile (110) surface. 6. 2 .2 Reconstructed Rutile (110) Surfaces In experiments, the rutile (110) surface is stable but it reconstructs at higher temperatures such as 500K to 800K. 151 The reconstructed rutile (110) surface has been the subject of considerable interest since many reconstructed geometries have been observed with experimental techniques such as low energy electron diffraction/Auger electron spectroscopy 175 and scanning tunneling microscopy 176, 177 In this section, the surface energies of three different reconstructed configurations of rutile (110) surfaces are calculated. Snapshots of these relaxed surfaces, along with the pristine (110) surface, are illustrated in Figs. 6 2 (a) to (d). The pristine (110) rutile TiO 2 surface, given in F igure 6 2 (a), is stable under normal temperatures and pressures. The first Figure 6 2 (b), in which one half of the bridging oxygen atoms on the surface are Ti 2 O 3 Figure 6 2 (c), in which an extra Ti 2 O 3 is added to the surface between model, which is provided in Figure 6 2 (d). In this case, a n extra row of Ti 3 O 5 is added to the surface. Surface energies of these reconstructed surfaces are calculated at two different charge states: neutral and charged. Since TiO 2 can be easily reduced, these reconstructed surfaces all are reduced relative to bulk TiO 2 ; neutral and charged surfaces therefore assume the missing oxygens are O ions and O atoms, respectively. Surface Ti atoms of the charged surfaces have approximately an additional 0.08 e of charge relative to the
130 neutral surfaces. The relaxed structures provided in Fig ure 6 2 are charged surfaces, and the surface energies of these reconstructed surfaces from COMB potential are compared to our DFT calculation values in Table 6 8 The DFT calculations were conducted with the same settings as described in Section 6. 1 .3 wit h the exception of a 441 Monkhorst Pack k point mesh is used. Snapshots of relaxed reconstructed surfaces with DFT are also illustrated in Figure 6 2 For un relaxed surfaces, both neutral and charged surface energies from COMB potential are predicted to be in the same order as that from DFT calculations. The pristine (110) surface is predicted to be more stable than any reconstructed (110) surface, consistent with the experimental finding that the pristine (110) surface is stable at lower temperatures. 151 For reconstructed surfaces, the missing row reconstruction is the most stable, followed by added Ti 2 O 3 and extra row reconstructions, indicating missing row reconstruction should be the most observed reconstructed surface. For COMB potential predictions, neutral reconstructed surfaces have lower surface energies than charged surfaces; this is because when O atoms, instead of O ions, are removed to form charged reconstructed surfaces, surface Ti atoms have excessive positive charges which adds extra energy to the overall potential energy from the contribution from the self energy term ( Section 2.2.1, Equatio n 2 1 2 ). This potential model artifact is physical since creating a charged reconstructed surface involves removing O atoms, instead of O ions, which requires more energy. Comparing COMB and DFT predictions, we find the over estimation in surface energies is consistent with the over estimation s predicted for defect formation energies.
131 For relaxed surfaces, on the other hand, we find the neutral reconstructed (110) surface energies are substantially lower as predicted by the COMB potential than the values f or charged surfaces For some cases the surfaces are even lower in energy than in bulk rutile. The unreasonably low relaxed surface energies result from the surface Ti atoms being significantly relaxed into the bulk position, in a similar fashion to the (1 00) surface discussed in the previous S ection ( 6. 2 .1 ) For charged surfaces, the excessive charges on surface Ti atoms not only increased the self energy ( Equation 2 1 2 ) but also increased the Coulomb repulsion between surface and bulk Ti atoms; hence the charge prevented the surface Ti atoms from being overly relaxed into the bulk. Comparing the relaxed structures from COMB potential (charged surfaces) and DFT, we find similarities between the two methods for the pristine surface, Figure 6 2 (a), the missi ng row reconstructed surface, Figure 6 2 (b), and the extra row reconstructed surface, Figure 6 2 (d). In addition, the trends for relaxed surface energies agree with each other, indicating that the pristine (110) surface is preferred compared to the missi ng row and the extra row reconstructions. However, the relaxed structures for added Ti 2 O 3 reconstruction are significantly different from DFT and COMB potential. From COMB potential, the 6 fold coordinated surface Ti atoms stay 6 fold, while from DFT those Ti atoms become 4 fold coordinated. The change in coordination number of the Ti 2 O 3 surface unit from DFT calculations results in the too low relaxed surface energy. In sum, from COMB potential and DFT calculations, we find the pristine (110) surface shoul d be more energetically preferred than the reconstructed (110) surfaces considered.
132 6. 2 .3 Step Rutile (110) Surfaces Rutile surfaces are often used as a heterogeneous catalyst which requires a certain amount of surface preparation. Repeated cycles of ultr a high vacuum (UHV) annealing and sputtering treatment is often used to clean surfaces. 178 After treating the surfaces under certain conditions (110) rutile can form a flat surface with (1x1) step terraces. At annealing temperatures of 1100 K the terrace width is typically 100 and will increase with annealing temperatures. The monatomic step height is approximately 3.2 179 The presence and geometries of steps in the rutile surface greatly affect the catalysis process wherein the deposited constituent will preferentially deposit near and on steps, 180 which is why i t warrants further investigation here. In this section, the surface energies were investigated for a total of four different step geometries. Visualization of the differing step surfaces can be seen in Figure 6 3 (a) (d). Surface steps on (1x1) rutile ru n predominately parallel to the  and directions. There are two distinct geometries among the steps in the  directions, one of which has a smooth termination along the step edge. The smooth  step can either be terminated by oxygen, Figu re 6 3 (a) or titanium, Figure 6 3 (b). Both geometries were tested here. The other distinct geometry along  is atomically rough and is thought to be a reconstruction of the smooth termination, shown in Figure 6 3 (c). It is observed that the rough an d smooth  steps occur at the same probability. 179 Lastly is the step surface, Figure 6 3 (d), which has been reported by STM studies that steps dominate in those direction. 179 All step surfaces were constructed by cutting away parts of the surface monolayer of a perfect (110) rutile surface, and were based on models proposed in Diebold [ 181 ].
133 Since the sizes of the step surfaces are too large to be considered with DFT calculations, the surface energies are calculated only with the COMB potential d eveloped in this work. The surface energies of these step (110) surfaces predicted with the COMB TiO 2 potential are presented in Table 6 9 In general, step configurations form on higher index surfaces in an attempt to revert to lower index surfaces. 178 Additionally, the COMB potential predicts that step configurations are also stable on the low index (110) rutile surface. As presented in Table 6 9 smooth step surfaces ha ve lower relaxed surface energies than the pristine (110) rutile surface. As the surface relaxes the titanium atom closest to the st While titanium relaxes into the surface, the oxygen titanium oxygen bond facing into the surface increases from approximately 60 o to 98 o which is closer to the ideal bond angle of 90 o This effectively reduces the bond angle repulsive energy (described in Equation 2 30) per O Ti O bond angle along the step edge. The reduction in this repulsive energy overpowers the normal increase in surface energy caused by step geometries due to under coordinated atoms and result s in a lower surface energy overall. There is still some controversy over the exact structure of the  smooth step on rutile (110) surfaces due to different interpretations of STM images. 181, 182 In this study the COMB potential predicts the oxygen terminated structu re should be the more stable model for the  smooth step surface due to its lowe r surface energy. The titanium terminated model has a higher surface energy because of increased number of under coordinated titanium atoms on the surface as compared to th e oxygen terminated  step surface. This is because t he surface energy is directly proportional to the
134 quantity of under coordinated atoms or the severity of under coordination, referring to the increase in excess energy per site as coordination decrea ses. 183 Although the rough  and smooth  step models appear approximately at the same probability experimentally 181, 182 t he rough  step has a higher surface energy predicted from the COMB TiO 2 potential as compared to the  smooth step models and the pristine (110) due to the under coordinated titanium atoms that protrude from the step Experimentally the step surface is more prevalent in STM studies at lower temperatures 181, 182 than the  steps and should have a lower surface energy. The higher surface energy predicted from COMB mentioned in Table 6 9 m ay be attributed to the interactions between the step edges. The surface area considered here may be too small hence allowing interference between the step edges resulting in a higher surface energy. In sum, COMB predicts that steps on rutile (110) surface s are more energetically favorable than the pristine (110) surface. 6. 2 4 Low Miller Index Anatase Surfaces Due to its surface properties, the anatase phase is believed to be a more efficient material than the rutile phase in applications such as catalysis photo catalysis and dye sensitized solar cells. 172 Therefore we assess the capability of COMB Ti/TiO 2 potential to describe anatase surfaces. Mentioned in Table 6 10 are the surface energies of the (101), (100), (001), (110), faceted (103) and smooth (103) surfaces of anatase as predicted from the COMB potential compared to first principles results reported by Lazzeri et al. 172, 173 Hummer et al. 171 Barnard et al. 170 and the MS Q 152 potential. Compared to values from Lazzeri et al, the un relaxed surface energies from COMB are typically over est imated by a factor of two for the similar reason to the rutile surface energies, which is the strong Coulombic contribution resulted from large fractional
135 charges. The MS Q potential correctly predicted the order between (100) and (001) surfaces, however n o information on the properties of anatase surfaces using SMB Q potential have been published to date. The relaxed surface energies from the COMB potential generally follow the DFT predicted trend. However the relaxed surface energies of (100) and smooth ( 103) surfaces are substantially higher compared then the other surfaces. The (101) surface is correctly predicted to be the most stable anatase surface. In conclusion, COMB Ti/TiO 2 potential generally captures the surface properties of TiO 2 polymorphs in t erms of surface energies and geometries. 6. 3 Adsorption Energies of Cu n (n=2~4) Clusters on Rutile (110) Surface In previous Sections 6. 1 and 6. 2 we find some empirical potentials out performs the COMB Ti/TiO 2 potential in some aspects, however due to thei r limitations in functional formalisms and transferabilities, it is difficult to model metal clusters interacting with titania surfaces. Therefore in this section we demonstrate the capability of COMB potentials to model such interactions. Because of the s tability of rutile titania (110) surface, it is able to support an extensive assortment of metal clusters including gold, lithium, aluminum, silver and copper. 151, 184 186 TiO 2 supported metal clusters are important catalysts for NO x reduction 187 189 and methanol production. 190, 191 Shown in Figure 6 4 is a top view of the (1x1) ru tile (110) surface indicating the adsorption sites. A first principles DFT study by Pillay et al. 192 compared the adsorption energies of Cu/Au dimers, trimers, an d tetramers with various geometries on specific adsorption sites on a rutile (110) surfaces. It is found for (1x1) surface on top of the Ti (5c) site is the most preferable adsorption site for Cu dimers/trimers and the structure of these Cu clusters is aff ected by the particle ionization by charge transfer from the TiO 2 substrate and Cu TiO 2 bonding interaction. Here we apply the COMB potentials for Ti/TiO 2 and Cu/Cu 2 O
136 to investigate the energetic and behavior of Cu n (n=2 4) clusters on the (1x1) rutile (11 0) surface. 6. 3 .1 Adsorption Energies and Geometries For the purpose of comparison, we considered the same rutile (110) supported Cu clusters as shown in Pillay et al. 192 for adsorption energy calculations. Force and energy are minimized followed by constant volume and constant temperature relaxation at 1K with charge equilibration. The relaxed geometries of adsorbed Cu dimers, trimers, and tetramers from COMB potential are illustrated in Figure 6 5 and the Cu Cu and Cu O bond length are summarized in Table 6 11 The adsorption geometries are generally comparable to those obtained with DFT calculations, 192 and all a dsorbed Cu clusters remain stable on the adsorption sites. Examining the relaxed adsorption geometries, we see that the Cu Cu bond lengths are typically underestimated by 15%, while the Cu O bond lengths are overestimated by ~25%. The adsorption energies a re also summarized in Table 6 11 The trend of these adsorption energies from COMB potentials generally follows that from DFT calculations: for Cu dimers the COMB potentials predict the same trend compared to DFT calculations. The Ti (5c) is the most prefe rred adsorption site while Cu dimer on top of O (2c) site with the bond perpendicular to  direction is the least preferred. For Cu trimers, the COMB potentials predict the same qualitative trend for Ti (5c) and O (2c) sites, however the O (3c) site is significantly underestimated so that this site is incorrectly predicted to be the least preferred adsorption site for Cu trimers. The adsorption energies of Cu tetramers from COMB potentials are also consistent with that of Cu dimers and trimers Ti (5c) is more preferred than O (3c) site. The results show that for Cu clusters the Ti (5c) is the most favored adsorption site, which is in agreement
137 with the finding of Pillay et al. and indicate COMB potentials should be adequate in describing the adsorptio n energetic of Cu clusters on rutile (1x1) (110) surface. The adsorption energies from COMB potentials are generally overestimated compared to that from DFT calculations with a largest overestimation of 60%. The overestimation results from the stronger Cu Ti short range attractions. 6.3. 2 Potential Energy Surface of a Cu Dimer on Rutile (110) Surface It is reported that Cu atoms/clusters have high mobility on rutile (110) surface at room temperature and even at temperatures as low as 160 K. 193, 194 Ther efore we applied the COMB potentials to calculate the potential energy surface of a Cu dimer moving along the  direction for a displacement of one surface unitcell on the (1x1) rutile (110) surface. The TiO 2 substrate is held fixed throughout the surf ace scan but is allowed to transfer charges with the Cu dimer. The charge equilibration is carried out with a 0.001 convergence criterion. The potential energy surface scan profile is shown in Figure 6 6 The potential energies (in J/mol) are plotted again st the lowest energy value, which corresponds to that of the Cu dimer placed on top of the Ti (5c) site. The minimum at the distance of 1.5 corresponds to the Cu dimer placed on top of the O (3c) site. The minima at the distances of ~0.5 (also 2.5 ) a nd 1.0 (also 2.0) correspond to the breaking of Cu Ti (5c) bond and Cu O (2c) bond, respectively, for each Cu atom. The breaking of these bonds decreased the coordination number of one Cu atom from 7 to 6, resulting in a slightly more preferred bonding environment for the Cu atom. The energy barrier for this Cu dimer moving on the rutile surface is estimated to be around 2.0 kJ /mol (~0.041 eV/dimer), and corresponds to a thermal energy of 480 K. This calculated energy barrier corresponds to the overesti mated adsorption energies and is reasonable compared to the experimental findings. 193, 194
138 6. 4 Summary To conclude, we have developed an empirical, variable charge potential for Ti and TiO 2 systems that is based on the charge optimized many body (COMB) potential framework. We also carried out first principles DFT calculations to determine the phase order of titania polymorphs and the elastic constants of brookite phase. The developed potential is fitted to, and captures most of, the structural, energetic and mechanical properties of Ti metal and various TiO 2 polymorphs. In addition, w e predicted defect formation energies and surface energies with the potential. The deficiencies of the potential include Ti surface energies, Ti interstitial formation energy certain anatase and brookite elastic constants and larger energy differences for several high pressure and hypothetical titania phases compared to DFT calculations. The potential generally reproduces the surface energies of several low Miller index rutil e, reconstructed rutile (110), and anatase surfaces. Combining the Ti/TiO 2 potential developed in this work with the Cu/Cu 2 O 30 potential, we found the adsorption energies of Cu dimers, trimers, and tetramers on rutile (110) surface are generally o verestimated compared to DFT calculations, however the adsorption geometries and trend of the adsorption energies are in reasonable agreement. The estimated energy barrier of a Cu dimer moving on the rutile (110) surface is also in agreement with experimen tal findings. Thus the COMB potential for Ti/TiO 2 should prove to be adequate in modeling bulk and surface properties of TiO 2 polymorphs. Additionally, the adsorption energetic and geometries of titania supported Cu clusters can be modeled in large scale a tomic level simulations.
139 Table 6 1. Properties of Ti metal given by COMB Ti/TiO 2 potential in comparison with experimental, DFT calculations and a MEAM potential Properties Exp DFT + MEAM + COMB a 0 () 2.951 2.948 2.931 2.933 c 0 () 4.684 4.667 4.678 4.754 c 0 /a 0 1.587 1.583 1.596 1.620 E c (eV/atom) 4.844 5.171 4.831 5.077 C 11 (GPa) 176 160 172 174 171 C 12 (GPa) 87 160 82 95 88 C 13 (GPa) 68 160 75 72 79 C 33 (GPa) 191 160 190 190 186 C 44 (GPa) 51 160 45 58 34 C 66 (GPa) 28 ++ 24 B (GPa) 110 111 113 113 G (GPa) 52 47 55 38 E ( fcc hcp ) (meV/atom) 60 159 58 39 58 E ( bcc hcp ) (meV/atom) 70 159 108 111 67 Vacancy (eV/defect) 2.03 2.24 1.79 Interstitial (eV/defect) 2.58 2.64 4.06 I 1 stacking fault (mJ/m 2 ) 124 I 2 stacking fault (mJ/m 2 ) 300 158 287 170 248 Extrinsic stacking fault (mJ/m 2 ) 227 Unstable stacking fault (mJ/m 2 ) 394 195 371 2 ) 1920 196 1939 1474 1260 ( mJ/m 2 ) 2451 1554 1308 (mJ/m 2 ) 1875 1682 1643 Denotes fitted properties for the COMB potential + DFT and MEAM values from Ref. [ 158 ] unless otherwise specified ++ This work
140 Table 6 2. Properties of rutile, anatase and brookite titania phases given by the COMB potential for Ti/TiO 2 developed in this work in comparison with that from experiments, our first principles calculations (values without superscripted references), a tight binding method, and three other empirical potentials Rutil e TiO 2 Properties Exp DFT COMB TB 157 SMB Q 1 56 MS Q 152 Matsui a 0 () 4.594 4.646 4.594 4.594 4.581 4.587 4.493 c 0 () 2.959 2.970 2.961 2.958 2.966 2.958 3.008 E c (eV/TiO 2 ) 26.95 26.77 19.91 19.90 16.49 39.80 B (GPa) 218 197 205 198 211 199 199 219 228 229 237 G (GPa) 124 197 107 198 136 199 125 113 87 116 C 11 (GPa) 268 197 276 198 270 199 299 290 293 294 322 C 12 (GPa) 175 197 154 198 172 199 156 178 203 202 230 C 13 ( GPa) 147 197 152 198 147 199 121 160 164 168 147 C 33 (GPa) 484 197 483 198 468 199 39 8 399 400 423 444 C 44 (GPa) 124 197 113 198 116 199 130 118 128 96 123 C 66 (GPa) 190 197 211 198 216 199 166 167 183 190 226 q Ti (e) 2.26 3.22 2.03 2.51 1.15 2.196 Anatase TiO 2 Properties Exp DFT COMB SMB Q 156 MS Q 152 Matsui a 0 () 3.785 3.806 3.785 3.825 3.85 3.770 c 0 () 9.512 9.714 9.514 9.03 9.06 9.568 E c (eV/TiO 2 ) 27.03 26.70 19.84 8.19 39.49 B (GPa) 59 178 152 180 199 164 220 176 176 G (GPa) 57 199 103 66 59 C 11 (GPa) 336 199 354 496 451 C 12 (GPa) 139 199 53 83 79 C 13 (GPa) 136 199 87 118 103 C 33 (GPa) 192 199 310 206 207 C 44 ( GPa) 49 199 95 51 25 C 66 (GPa) 58 199 60 54 49 q Ti (e) 2.24 3.23 1.12 2.196 Brookite TiO 2 Properties Exp DFT COMB SMB Q 156 MS Q 152 Matsui a 0 () 9.184 9.267 9.160 9.259 9.113 9.146 b 0 () 5.449 5.505 5.430 5.444 5.449 5.389 c 0 () 5.145 5.180 5.130 5.229 5.170 5.144 E c (eV/TiO 2 ) 26.99 26.62 19.65 8.21 39.62 B (GPa) 179 152 192 177 227 211 200 G (GPa) 125 102 83 92 C 11 (GPa) 274 264 272 262 C 12 (GPa) 141 132 177 183
141 Table 6 2. Continued Properties Exp DFT COMB SMB Q 156 MS Q 152 Matsui 155 C 13 (GPa) 165 114 165 165 C 22 (GPa) 278 288 319 332 C 23 (GPa) 136 97 164 124 C 33 (GPa) 291 348 301 273 C 44 (GPa) 101 97 90 99 C 55 (GPa) 104 114 97 113 C 66 (GPa) 90 177 141 158 q Ti (e) 2.25 3.23 1.14 2.196 Denotes fitted properties for the COMB potential
142 Table 6 3. Cohesive energies of rutile titania (eV / TiO 2 unit) and that relative to rutile titania for other polymorphs (including hypothetical AB 2 polymorphs) calculated from COMB potential for Ti / TiO 2 developed in this work in comparison with our first principles DFT calculations and other em pirical potentials Polymorph (space group) DFT COMB SMB Q 156 MS Q 152 Matsui 155 Normal phases Rutile (P42/mnm) 26.949 26.778 19.90 8.248 39.798 Anatase (I41/amd) 0.083 0.075 0.060 0.057 0.304 Brookite (Pbca) 0.046 0.162 0.250 0.037 0.178 High pressure phases Columbite (Pbcn) 0.006 2.564 0.074 0.035 Cotunnite (Pnma) 0.107 4.476 0.553 0.914 Baddeleyite (P21/c) 0.066 0.844 0.515 Fluorite (Fm 3m) 0.775 1.031 1.020 0.580 Hypothetical phases Quartz (P3221) 0.228 0.798 Quartz (P6222) 4.296 1.512 Cristobalite (P41212) 0.108 0.587 Cristobalite (Fd 3m) 0.106 0.851 Marcasite (Pnnm) 0.016 2.679 HgCl2 (Pnma) 0.035 2.989 Table 6 4. Defect formation energies of bulk rutile Defect Defect formation energies (eV/defect) COMB Potential Exp. 163 DFT 164 Matsui Oxidized ( + ) Standard ( + ) Reduced ( + ) Ti interstitial 7.89 1.57 1.57 2.6 0.91 19.90 O vacancy 2.09 2.09 3.02 2.1 1.19 9.45 O interstitial 15.89 15.89 21.00 -3.06 0.14 Ti vacancy 28.64 34.97 34.97 2.4 4.94 38.44 Schottky 29.99 36.32 26.10 6.6 8.50 -Cation Frenkel 33.02 5.0 7.05 -Anion Frenkel 15.88 -6.90 -
143 Table 6 5. Defect formation energies at rutile grain boundaries. Defect Defect formation energies (eV/defect) COMB Potential DFT 169 Oxidized Standard Reduced Schottky (mixed A ) 14.88 25.26 4.03 3.11 Schottky (clustered) 29.27 39.65 18.42 3.19 Schottky (mixed B ) 30.28 40.66 19.43 5.72 Schottky (distributed) 30.00 40.38 19.15 6.02 Cation Frenkel 6.33 3.05 Anion Frenkel 13.37 12.30 Table 6 6 Coordination number analysis of Ti and O atoms from the rutile structures numbers are percentages of the Ti or O atom (the coordination number indicated in parenthesis) relative to total Ti or O atoms, respectively. The coordinated (6 for Ti, 3 for O). Defects O (2) O (3) O (4) Ti (4) Ti (5) Ti (6) Normal Schottky (mixed A ) 14.3 83.3 2.4 0.0 23.8 76.2 81.0 Schottky (clustered) 14.3 85.7 0.0 4.8 19.0 76.2 82.5 Schottky (mixed B ) 19.0 78.6 2.4 7.1 19.0 73.8 77.0 Schottky (distributed) 16.7 83.3 0.0 2.4 28.6 69.0 78.6 Cation Frenkel 9.1 88.6 2.3 0.0 18.2 81.8 86.4 Anion Frenkel 11.4 88.6 0.0 0.0 27.3 72.7 83.3 Table 6 7 Surface energies of rutile surfaces (J/m 2 ) COMB Ramamoorthy 174 Humme r 171 Barnard 170 Lazzeri 172, 173 MS Q 152 SMB Q 156 Un relaxed Relaxe d Un relaxed Relaxed Relaxed Relaxe d Relaxe d Relaxed Relaxe d (110) 0.94 0.67 1.76 0.90 0.69 0.47 0.31 -0.68 (100) 1.34 0.53 1.94 1.12 ---2.17 0.85 (101) 1.72 0.98 2.12 1.40 1.12 0.95 ---(001) 3.24 1.55 2.95 1.65 ---1.86 1.74 Table 6 8 Surface energies of various reconstructed rutile (110) surfaces Rutile (110) surfaces COMB (neutral) COMB (charged) DFT Un relaxed Relaxed Un relaxed Relaxed Un relaxed Relaxed Pristine 0.94 0.67 0.94 0.67 0.97 0.50 Missing row 1.30 0.54 2.22 1.14 1.26 0.51 Added Ti 2 O 3 4.42 0.05 6.18 1.35 5.10 0.15 Extra row 9.15 0.20 9.51 2.44 5.82 0.55
144 Table 6 9 Surface energies of rutile step (110) surfaces Step surface Surface energy from COMB (J/m 2 ) Smooth  O 0.4 7 Smooth  Ti 0.5 8 Rough  0.9 5 1. 12 Pristine (110) 0.67 Table 6 10 Surface energies of anatase surfaces (J/m 2 ) COMB Lazzeri 172, 173 Hummer 171 Barnard 170 MS Q 152 Un relaxed Relaxed Un relaxed Relaxed Relaxed Relaxed Relaxed (101) 2.44 0.59 1.28 0.44 0.41 0.35 -(100) 3.15 1.18 1.59 0.53 0.51 0.39 1.34 (103) f 3.79 0.71 1.50 0.83 ---(001) 4.01 0.73 1.12 0.90 0.96 0.51 2.04 (103) s 5.60 1.72 2.40 0.93 ---(110) 4.46 1.34 2.17 1.09 ---
145 Table 6 11 Adsorption bond lengths and energies (in bold) of Cu dimers, trimers, and tetramers on the (1x1) rutile TiO 2 (110) surface. A B C D Cu 2 DFT 192 Cu Cu () 2.42 3.02 2.36 2.53 Cu O () 2.04 1.88 1.86 2.02 E ad (eV) 3.07 2.57 2.54 2.27 COMB Cu Cu () 2.19 3.23 2.19 2.09 Cu O () 2.54 1.90 2.54 2.55 E ad (eV) 4.30 4.26 4.20 3.28 Cu 3 DFT 192 Cu Cu () 2.32, 2.66 2.33, 2. 61 2.36, 2.68 2.94 Cu O () 2.00 2.46 2.10 2.10 E ad (eV) 3.61 3.30 3.11 1.63 COMB Cu Cu () 2.14, 2.27 2.15, 2.27 2.14, 2.27 3.00 Cu O () 2.55 2.37 2.15 2.15 E ad (eV) 4.89 1.36 3.31 2.15 Cu 4 DFT 192 Cu Cu () 2.37 2.36~2.43 -Cu O () 2.38 2.50 E ad (eV) 3.69 3.30 COMB Cu Cu () 2.21 2.23~2.39 Cu O () 2.66 2.65 E ad (eV) 4.19 1.60
146 A B C Figure 6 1. grain boundaries relaxed with the COMB potential developed in this work: A ) Schottky (mixed A configuration), B ) Schottky (clustered configuration), C ) Schottky (mixed B configuration ), D ) Schottky (distributed configuration), E ) Cation Frenkel, and F ) Anion Frenkel defects. Larger atoms are Ti; smaller atoms are O. Atoms are color coded by the coordination number: maroon=6, red=5, gray=4, green=3 and wood=2.
147 D E F Figure 6 1. Continued
148 A B Figure 6 2. Pristine and reconstructed rutile TiO 2 (110) surfaces relaxed with the COMB potential developed in this work (left) and first principles DFT calculations (right): A ) pristine, B ) missing row reconstructed, C ) adde d Ti 2 O 3 reconstructed, and D ) extra row reconstructed surfaces. Larger atoms are Ti; smaller atoms are O. Atoms are color coded by the coordination number: maroon=6, red=5, gray=4, green=3 and wood=2.
149 C D Figure 6 2. Continued
150 Figure 6 3. Relaxed geometries of rutile (110) step surfaces with the COMB Ti/TiO 2 potential developed in this work. A ) Smooth O terminated , B ) smooth Ti ter minated , C ) rough , and D ) step surfaces. Large gray spheres are Ti, small red spheres are O; surface Ti atoms are highlighted large yellow spheres.
151 Figure 6 4 Sphere representation of the (1x1) rutile TiO 2 (110) surface viewing from the surface normal direction. Larger atoms are Ti; smaller atoms are O. Atoms are color coded by the coordination number: maroon=6, red=5, green=3 and wood=2. Ti (6c) Ti (5c) O (3c) O (2c)
152 A B C Figure 6 5 Adsorption geometries of Cu A ) dimers, B) trimers and C) tetramers on the (1x1) rutile TiO 2 (110) surface. Atoms are color coded by the charge values indicated by the color bar. Ti(5c) O(3c) Ti(5c) O (3c) O(2c) O(2c) Ti(5c) O(2c) Ti(6c) O(2c)
153 Figure 6 6 Potential energy surface of a Cu dimer moving on (1x1) rutile (110) surface along  direction for one surface unitcell. O(3c) Ti(5c) O(2c) Ti(5c)
154 CHAPTER 7 CONCLUSIONS 7.1 General Summaries An empirical, variable charge many body (COMB) potential has been developed for modeling multifunctional nanostructures composed of heterogeneous interfaces. The first generation COMB potential function was designed to model the discrete bonding of covalent, metallic and ionic bonds in an integrated scheme. 12, 40 In this work an improved, second generation COM B potential is developed for Si/SiO 2 16 Cu/Cu 2 O 37 Hf/HfO 2 13 and Ti/TiO 2 38 systems. The achievements of the second generation COMB Si/SiO 2 potential include the replacement of the unphysical treatment of long range Coulombics and improved the mechanical properties and phase stabil ity of oxide polymorphs over the first generation COMB potential. Another main achievement of the second generation COMB potential is its use to describe amorphous oxide phases, which the first generation COMB potential was unable to do. The oxygen paramet ers are exactly the same for SiO 2 Cu 2 O, HfO 2 and TiO 2 potentials. This is an important feature since it allows all elemental and oxide phases to be modeled in the same simulation cell. Therefore it is possible to model complex multifunctional nanostructur es composed of heterogeneous interfaces in an integrated manner. These second generation COMB potentials can be used in molecular dynamics simulations to investigate the structural, energetic and mechanical properties of the fitted semiconductor, metal and metal oxide materials, as well as the adhesive properties of their interfaces. In this work successful application of the COMB potentials to model heterogeneous interfacial systems include:
155 1. Atomic oxygen adsorption energetic on reconstructed Si (001) surfaces, 2. Si nanocrystals embedded in amorphous silica matrix, 3. Nanoindentation simulations of the Si/ quartz, Si/a SiO 2 and Si/a HfO 2 systems, 4. Adhesion of Cu/ cristobalite, Cu/ quartz and Cu/a SiO 2 interfaces, and 5. Adsorption of Cu dimers, trimers and tetramers on rutile TiO 2 (110) surface. The empirical, variable charge second generation COMB potentials developed in this work are useful and powerful tools that will allow scientists and engineers to play in as described in Figure 1 1. With continued advancement in computer technology, real device size multifunctional nanostructures can be modeled at the atomistic level with MD simulations using the COMB potentials, thus providing insights into the atomic le vel properties of materials with technological significance. 7 2 Ongoing Development of COMB Potentials An important aspect of an empirical potential is its computational efficiency. As discussed in Section 2.5 and presented in Table 2 1, the second genera tion COMB potential is a relatively computational time intensive empirical potential function compared to others. For the COMB potential, approximately 75 % of the computational time is taken by charge equilibration. Provided that charge equilibration effic iency is improved, for example by 30%, the overall computational efficiency of the COMB potential can be increased by 25%. Currently another charge equilibration method, the iterative matrix method, is being investigated. The objective of this alternative method is to increase the charge equilibration efficiency while at least retain the precision. Another approach to increase the overall computational efficiency is to utilize look up tables for the values for the bond order term. Instead of calculating the bond order term (Section 2.2, Equation 2 18) at each time step, the bond order term is only calculated at the first
156 time step and a look up table for subsequent time steps is generated. From preliminary tests by Devine, the short range charge independent calculation efficiency is increased by a factor of three. This method is currently being investigated to be implemented in LAMMPS software. The second generation COMB potentials provide insights and foundations for the development of new potential function s for new elements and compounds. The author is participating in the development of COMB potentials for U/UO 2 Zr/ZrO 2 Zn/ZnO and Ti/TiN/TiO 2 Other new COMB potentials include potentials for C/H/O systems. With the additions of these new potential funct ions, a wide range of multifunctional nanostructures and heterogeneous interfaces of significance and importance can be investigated with the COMB potentials.
157 APPENDIX POTENTIAL PARAMETERS FOR COMB POTENTIALS DEVELOPED IN THIS WO RK The potential parameters for COMB potentials for Si/SiO 2 Cu/Cu 2 O, Hf/HfO 2 and Ti/TiO 2 developed in this are given in Table A 1. The interaction cutoff radii and mixing rule scaling factors are given in Table A 2.
158 Table A 1. Potential parameters for COMB potentials developed in this work Si O Cu Hf Ti A (eV) 1803.80 3326.699 952.693 707.5303 546.3860 B (eV) 471.18 260.8931 146.987 55.94216 99.3916 ( 1 ) 2.4799 5.36 2.794608 2.069563 2.0824 1 ) 1.7322 2.68 1.681711 0.959614 1.2263 1.0999 10 6 2.0 0.140835 0.046511 0.089078 n 0.78734 1 1 1.011011 1.255016 m 3 1 1 1 1 c 100390 6.6 0 0 0 d 16.218 1 1 1 1 h 0.59826 0.229 0 0 0 R S () 2.8 2.6 3.1 5 3. 8 3.27 S S () 3.0 3.0 3.3 5 4.0 3.57 Q L 4.0 1.8349 6.0 ( 7.2) 4.0 4 Q U 4.0 5.5046 2.0 (0.8) 4.0 4 D L 1.651725 0.00148 0.16776 0.26152 2.508854 D U 1.658949 0.00112 0.16100 0.25918 2.511416 n B 10 10 10 10 10 0 5.63441 0 0 0 J(eV) 3.62514 7.68960 5.946437 3.13952 2.46820 K(eV) 0 4.51427 0 0 0 L(eV) 0.08707 1.33008 0 0.00941 0.15135 0.772871 2.243072 0.454784 0.679131 0.873685 1 0.499378 3.922011 0.72571 3.928750 0.392632 2 2.999911 0.971086 0.274649 4.839580 1.78349 (eV) 2.60 -0.007858 -0.202777 (eV) 0.3122 -2.518789 -0.403105 143.73 -109.47 -130.54 109.47 -180.0 -90.0 --0.077 ----0.0095 -----0.008 0.0084 K r ----8.45 E 0 ---0.16 ----0.10 -R ---0.14 -
159 Table A 2 Cutoff radii and m ixing rule scaling factors for short range interactions. Si O Si Cu Cu O Hf X T i Si T i O Ti Cu R () 2.55 2.80 2.45 3.17 -3.15 4.35 S () 3.05 2.90 2.65 3.41 -3.35 4.55 A 1.0 1.276957 1.666072 1.0 1.0 1.326746 10.000363 B 1.0 0.452693 0.100000 1.0 1.0 0.088406 0.0981790 1.0 1.032042 1.097775 1.0 1.0 0.969934 1.6622380 1.0 1.354486 0.584713 1.0 1.0 0.296577 0.1048210
160 REFERENCE LIST 1 S. R. Phillpot, and S. B. Sinnott, Science 325 1634 (2009). 2 R. McWeeny, Methods of Molecular Quantum Mechanics (Academic, New York City, New York, 1989). 3 P. Hohenberg and W. Kohn, Phys. Rev. B 136 B864 (1964). 4 W. Kohn, and L. J. Sham, Phys. Rev. 140 1133 (1965). 5 M. S. Daw, and M. I. Baskes, Phys. Rev. B 29 6443 (1984). 6 S. M. Foiles, M. I. Baskes, and M. S. Daw, Phys. Rev. B 33 7983 (1986). 7 R. A. Buckingham Proceedings of the Royal Society of London, Series A: Mathematical and Physical Sciences 168 264 (1938). 8 J. Tersoff, Physical Review B 37 6991 (1988). 9 J. Tersoff, Physical Review B 38 9902 (1988). 10 A. C. T. van Duin, S. Dasgupta, F. Lorant, and W. A. Goddard, Journal of Physical Chemistry A 105 9396 (2001). 11 A. C. T. van Duin, A. Strachan, S. Stewman, Q. S. Zhang, X. Xu, and W. A. Goddard, Journal of Physical Chemistry A 107 3803 (2003). 12 J. Yu, S. B. Sinnott, and S. R. Phillpot, Physical R eview B 75 085311 (2007). 13 T. R. Shan, B. D. Devine, T. W. Kemper, S. B. Sinnott, and S. R. Phillpot, Physical Review B 81 125328 (2010). 14 M. P. Allen, and D. J. Tildesley, Computer Simulation of Liquids (Oxford University Press, New York, 1989). 15 J. E. Lennard Jones, Proceedings of the Royal Society of London. Series A 106 463 (1924). 16 T. R. Shan, B. D. Devine, J. M. Hawkins, A. Asthagiri, S. R. Phillpot, and S. B. Sinnott, Physical Review B 82 235302 (2010).
161 17 B. D. Devine, PhD Dissertation, University of Florida, 2010. 18 L. Verlet, Physical Review 159 98 (1967). 19 H. C. Andersen, Journal of Chemical Physics 72 2384 (1980). 20 W. G. Hoover, Physical Review A 31 1695 (1985). 21 S. Nose, Molecular Physics 52 255 (1984). 22 S. Nose, Journal of Chemical Physics 81 511 (1984). 23 H. J. C. Berendsen, J. P. M. Postma, W. F. Vangunsteren, A. Dinola, and J. R. Haak, J. Chem. Phys. 81 3684 (1984). 24 G. J. Martyna, D. J. Tobias, and M. L. Klein, J. Chem. Phys. 101 4177 (1994). 25 A. Yasukawa, Js me International Journal Series a Mechanics and Material Engineering 39 313 (1996). 26 A. K. Rappe, and W. A. Goddard, Journal of Physical Chemistry 95 3358 (1991). 27 F. H. Streitz, and J. W. Mintmire, Phys. Rev. B 50 11996 (1994). 28 T. Iwasaki, Journ al of Materials Research 19 1197 (2004). 29 T. Iwasaki, J. Mater. Res. 20 1300 (2005). 30 B. D. Devine, T. R. Shan, M. Y. Lee, A. McGaughey, S. B. Sinnott, and S. R. Phillpot, (submitted for review). 31 J. C. Slater, Physical Review 36 57 (1930). 32 D. Wolf, P. Keblinski, S. R. Phillpot, and J. Eggebrecht, Journal of Chemical Physics 110 8254 (1999). 33 D. W. Brenner, Physical Review B 42 9458 (1990). 34 D. W. Brenner, O. A. Shenderova, J. A. Harrison, S. J. Stuart, B. Ni, and S. B. Sinnott, Journal of Physics: Condensed Matter 14 783 (2002).
162 35 D. W. Brenner, Physica Status Solidi B: Basic Research 217 23 (2000). 36 M. W. Finnis, and J. E. Sinclair, Philosophical Magazine a Physics of Condensed Matter Structure Defects and Mechanical Properties 50 45 (1984). 37 T. R. Shan, B. D. Devine, S. R. Phillpot, and S. B. Sinnott, Physical Review B 83 115327 (2011). 38 T. R. Shan, R. K. Behera, J. A. Martinez, S. R. Phillpot, and S. B. Sinnott, (submitted for review). 39 B. J. Thijsse, Nuclear Instruments & Methods in Physics Research Section B Beam Interactions with Materials and Atoms 228 198 (2005). 40 J. Yu, S. B. Sinnott, and S. R. Phillpot, Philosophical Magazine Letters 89 136 (2009). 41 W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Fla nnery, Numerical Recipes in Fortran 77: The Art of Scientific Computing (Cambridge University Press, 1992). 42 Y. T. Cheng, B. D. Devine, T. R. Shan, S. R. Phillpot, and S. B. Sinnott, (in preparation). 43 Y. Li, T. R. Shan, S. B. Sinnott, and S. R. Phillp ot, (in preparation). 44 T. Liang, B. D. Devine, T. W. Kemper, T. R. Shan, S. B. Sinnott, and S. R. Phillpot, (in preparation). 45 M. J. Noordhoek, T. R. Shan, S. B. Sinnott, and S. R. Phillpot, (in preparation). 46 R. T. Sanderson, Science 114 670 (1951) 47 Z. X. Zhou, and R. G. Parr, Journal of the American Chemical Society 111 7371 (1989). 48 R. G. Parr, and W. Yang, Density Functional Theory of Atoms and Molecules (Oxford University Press, USA, 1994). 49 S. W. Rick, S. J. Stuart, and B. J. Berne, Jou rnal of Chemical Physics 101 6141 (1994).
163 50 S. J. Plimpton, Journal of Computational Physics 117 1 (1995). 51 F. H. Stillinger, and T. A. Weber, Physical Review B 31 5262 (1985). 52 E. C. T. Chao, E. M. Shoemaker, and B. M. Madsen, Science 132 220 (1960). 53 E. C. T. Chao, J. J. Fahey, J. Littler, and D. J. Milton, Journal of Geophysical Research 67 419 (1962). 54 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77 3865 (1996). 55 S. Tsuneyuki, M. Tsukada, H. Aoki, and Y. Matsui, Phys. R ev. Lett. 61 869 (1988). 56 B. W. H. van Beest, G. J. Kramer, and R. A. van Santen, Physical Review Letters 64 1955 (1990). 57 J. L. Rosenholtz, and D. T. Smith, Am. Mineral. 26 103 (1941). 58 R. L. Mozzi, and B. E. Warren, J. Appl. Crystallogr. 2 164 (1969). 59 P. G. Coombs, J. F. Denatale, P. J. Hood, D. K. McElfresh, R. S. Wortman, and J. F. Shackelford, Philosophical Magazine B Physics of Condensed Matter Statistical Mechanics Electronic Optical and Magnetic Properties 51 L39 (1985). 60 L. Martin S amos, Y. Limoge, J. P. Crocombette, G. Roma, N. Richard, E. Anglada, and E. Artacho, Phys. Rev. B 71 014116 (2005). 61 N. Richard, L. Martin Samos, G. Roma, Y. Limoge, and J. P. Crocombette, J. Non Cryst. Solids 351 1825 (2005). 62 G. Kresse, and J. Furt hmuÂ¨ller, Computational Materials Science 6 15 (1996). 63 G. Kresse, and D. Joubert, Phys. Rev. B 59 1758 (1999). 64 P. E. Blochl, Phys. Rev. B 50 17953 (1994). 65 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). 66 M. Methfessel, and A. T. Paxton, Phys. Rev. B 40 3616 (1989).
164 67 H. J. Monkhorst, and J. D. Pack, Physical Review B 13 5188 (1976). 68 Y. Widjaja, and C. B. Musgrave, J. Chem. Phys. 116 5774 (2002). 69 A. L Gerrard, J. J. Chen, and J. F. Weaver, J. Phys. Chem. B 109 8017 (2005). 70 R. F. W. Bader, W. H. Henneker, and P. E. Cade, J. Chem. Phys. 46 3341 (1967). 71 V. Vinciguerra, G. Franzo, F. Priolo, F. Iacona, and C. Spinella, Journal of Applied Physics 8 7 8165 (2000). 72 T. Takagahara, and K. Takeda, Physical Review B 46 15578 (1992). 73 L. W. Wang, and A. Zunger, Journal of Physical Chemistry 98 2158 (1994). 74 M. V. Wolkin, J. Jorne, P. M. Fauchet, G. Allan, and C. Delerue, Physical Review Letters 82 197 (1999). 75 L. Dal Negro, M. Cazzanelli, N. Daldosso, Z. Gaburro, L. Pavesi, F. Priolo, D. Pacifici, G. Franzo, and F. Iacona, Physica E Low Dimensional Systems & Nanostructures 16 297 (2003). 76 N. Daldosso, M. Luppi, et al., Physical Review B 68 (2 003). 77 F. Djurabekova, M. Backholm, M. Backman, O. H. Pakarinen, J. Keinonen, K. Nordlund, T. R. Shan, B. D. Devine, and S. B. Sinnott, Nuclear Instruments & Methods in Physics Research, Section B: Beam Interactions with Materials and Atoms 268 3095 (20 10). 78 B. R. Lawn, and D. B. Marshall, Journal of the American Ceramic Society 62 347 (1979). 79 J. B. Pethica, R. Hutchings, and W. C. Oliver, Philosophical Magazine A: Physics of Condensed Matter: Defects and Mechanical Properties 48 593 (1983). 80 W. C. Oliver, and G. M. Pharr, Journal of Materials Research 7 1564 (1992). 81 G. M. Pharr, W. C. Oliver, and F. R. Brotzen, Journal of Materials Research 7 613 (1992). 82 T. Chudoba, and F. Richter, Surface & Coatings Technology 148 191 (2001).
165 83 W. C. Oliver, and G. M. Pharr, Journal of Materials Research 19 3 (2004). 84 K. R. Morasch, and D. F. Bahr, Thin Solid Films 515 3298 (2007). 85 L. Kogut, and K. Komvopoulos, Journal of Materials Research 19 3641 (2004). 86 C. L. Kelchner S. J. Plimpton, and J. C. Hamilton, Physical Review B 58 11085 (1998). 87 V. B. Shenoy, R. Miller, E. B. Tadmor, D. Rodney, R. Phillips, and M. Ortiz, Journal of the Mechanics and Physics of Solids 47 611 (1999). 88 E. B. Tadmor, R. Miller, R. Phillips and M. Ortiz, Journal of Materials Research 14 2233 (1999). 89 A. Gouldstone, K. J. Van Vliet, and S. Suresh, Nature 411 656 (2001). 90 J. A. Zimmerman, C. L. Kelchner, P. A. Klein, J. C. Hamilton, and S. M. Foiles, Physical Review Letters 87 art. no. (2001). 91 J. Li, K. J. Van Vliet, T. Zhu, S. Yip, and S. Suresh, Nature 418 307 (2002). 92 E. T. Lilleodden, J. A. Zimmerman, S. M. Foiles, and W. D. Nix, Journal of the Mechanics and Physics of Solids 51 901 (2003). 93 H. P. Chen, R. K. Kalia, A. Naka no, P. Vashishta, and I. Szlufarska, Journal of Applied Physics 102 (2007). 94 I. Szlufarska, R. K. Kalia, A. Nakano, and P. Vashishta, Applied Physics Letters 85 378 (2004). 95 I. Szlufarska, R. K. Kalia, A. Nakano, and P. Vashishta, Applied Physics Lett ers 86 (2005). 96 I. Szlufarska, R. K. Kalia, A. Nakano, and P. Vashishta, Physical Review B 71 (2005). 97 I. Szlufarska, R. K. Kalia, A. Nakano, and P. Vashishta, Journal of Applied Physics 102 (2007). 98 I. Szlufarska, A. Nakano, and P. Vashishta Science 309 911 (2005).
166 99 P. Horvath, S. B. Sadale, et al., Sensor Letters 6 558 (2008). 100 P. Vashishta, R. K. Kalia, and A. Nakano, Journal of Physical Chemistry B 110 3727 (2006). 101 O. Yeheskel, R. Chaim, Z. J. Shen, and M. Nygren, Journal of M aterials Research 20 719 (2005). 102 C. W. Yong, K. Kendall, and W. Smith, Philosophical Transactions of the Royal Society of London Series a Mathematical Physical and Engineering Sciences 362 1915 (2004). 103 H. Kimizuka, H. Kaburaki, and Y. Kogure, Phy s. Rev. Lett. 84 5548 (2000). 104 T. Demuth, Y. Jeanvoine, J. Hafner, and J. G. Angyan, J. Phys.: Condens. Matter 11 3833 (1999). 105 S. P. Murarka, Mater. Sci. Eng., R 19 87 (1997). 106 M. Hasan, and J. F. Rohan, J. Electrochem. Soc. 157 D278 (2010). 107 S. P. Murarka, R. J. Gutmann, A. E. Kaloyeros, and W. A. Lanford, Thin Solid Films 236 257 (1993). 108 K. Nagao, J. B. Neaton, and N. W. Ashcroft, Physical Review B 68 125403 (2003). 109 J. Cechal, J. Polcak, M. Kolibal, P. Babor, and T. Sikola, Appl Surf. Sci. 256 3636. 110 Y. L. Cheng, T. J. Chiu, B. J. Wei, H. J. Wang, J. Wu, and Y. L. Wang, J. Vac. Sci. Technol., B 28 567 (2010). 111 M. D. Kriese, N. R. Moody, and W. W. Gerberich, Acta Mater. 46 6623 (1998). 112 M. Z. Pang, and S. P. Baker, J. Mater. Res. 20 2420 (2005). 113 T. S. Oh, R. M. Cannon, and R. O. Ritchie, Journal of the American Ceramic Society 70 C352 (1987). 114 J. Hallberg, and R. C. Hanson, Physica Status Solidi 42 305 (1970).
167 115 G. M. Rignanese, J. Phys.: Condens. Matter 17 R357 (2005). 116 X. Y. Zhao, and D. Vanderbilt, Phys. Rev. B 65 233106 (2002). 117 E. P. Gusev, C. Cabral, M. Copel, C. D'Emic, and M. Gribelyuk, Microelectron. Eng. 69 145 (2003). 118 H. Watanabe, S. Kamiyama, et al., Jpn. J. Appl. Phys., Part 2 44 L 1333 (2005). 119 N. Umezawa, K. Shiraishi, et al., Appl. Phys. Lett. 86 143507 (2005). 120 V. Cuny, and N. Richard, J. Appl. Phys. 104 (2008). 121 K. Tse, and J. Robertson, Appl. Phys. Lett. 89 (2006). 122 K. Xiong, and J. Robertson, Microelectron. Eng. 8 0 408 (2005). 123 P. F. Lee, J. Y. Dai, K. H. Wong, H. L. W. Chan, and C. L. Choy, J. Appl. Phys. 93 3665 (2003). 124 R. Choi, S. C. Song, C. D. Young, G. Bersuker, and B. H. Lee, Appl. Phys. Lett. 87 (2005). 125 R. M. Nieminen, M. H. Hakala, and A. S. F oster, Mater. Sci. Semicond. Process. 9 928 (2006). 126 R. C. Weast ed., Handbook of Chemistry and Physics (The Chemical Rubber Co., Cleveland, Ohio, 1969). 127 J. M. Leger, A. Atouf, P. E. Tomaszewski, and A. S. Pereira, Phys. Rev. B 48 93 (1993). 128 M Balog, M. Schieber, M. Michman, and S. Patai, Thin Solid Films 41 247 (1977). 129 P. K. Schelling, S. R. Phillpot, and D. Wolf, J. Am. Ceram. Soc. 84 1609 (2001). 130 M. Wilson, U. Schonberger, and M. W. Finnis, Phys. Rev. B 54 9147 (1996). 131 S. Fab ris, A. T. Paxton, and M. W. Finnis, Phys. Rev. B 61 6617 (2000).
168 132 X. Luo, A. A. Demkov, D. Triyoso, P. Fejes, R. Gregory, and S. Zollner, Phys. Rev. B 78 (2008). 133 J. E. Jaffe, R. A. Bachorz, and M. Gutowski, Phys. Rev. B 72 144107 (2005). 134 J. W ang, H. P. Li, and R. Stevens, J. Mater. Sci. 27 5397 (1992). 135 A. Christensen, and E. A. Carter, Phys. Rev. B 62 16968 (2000). 136 A. S. Foster, F. L. Gejo, A. L. Shluger, and R. M. Nieminen, Phys. Rev. B 65 174117 (2002). 137 M. Ritala, M. Leskela, L. Niinisto, T. Prohaska, G. Friedbacher, and M. Grasserbauer, Thin Solid Films 250 72 (1994). 138 S. V. Ushakov, A. Navrotsky, et al., Physica Status Solidi B 241 2268 (2004). 139 Q. Fang, J. Y. Zhang, et al., Thin Solid Films 453 203 (2004). 140 M. Y. Ho, H. Gong, et al., J. Appl. Phys. 93 1477 (2003). 141 D. Triyoso, R. Liu, et al., J. Electrochem. Soc. 151 F220 (2004). 142 H. Wong, and H. Iwai, Microelectron. Eng. 83 1867 (2006). 143 K. Cherkaoui, S. Monaghan, et al., J. Appl. Phys. 104 064113 (2008). 144 T. T. Tan, Z. T. Liu, H. C. Lu, W. T. Liu, and H. Tian, Optical Materials 32 432 (2010). 145 D. M. Adams, S. Leonard, D. R. Russell, and R. J. Cernik, J. Phys. Chem. Solids 52 1181 (1991). 146 D. W. Stacy, D. R. Wilder, and Johnston.Jk, J. Am Ceram. Soc. 55 482 (1972). 147 C. K. Lee, E. N. Cho, H. S. Lee, C. S. Hwang, and S. W. Han, Phys. Rev. B 78 012102 (2008). 148 J. E. Lowther, J. K. Dewhurst, J. M. Leger, and J. Haines, Phys. Rev. B 60 14485 (1999).
169 149 S. Desgreniers, and K. Lagarec, Phys. Rev. B 59 8467 (1999). 150 C. A. Ponce, R. A. Casali, and M. A. Caravaca, J. Phys.: Condens. Matter 20 045213 (2008). 151 U. Diebold, Surface Science Reports 48 53 (2003). 152 V. Swamy, J. D. Gale, and L. S. Dubrovinsky, Journal of Physics and Ch emistry of Solids 62 887 (2001). 153 D. R. Collins, W. Smith, N. M. Harrison, and T. R. Forester, Journal of Materials Chemistry 6 1385 (1996). 154 G. V. Lewis, and C. R. A. Catlow, Journal of Physics C Solid State Physics 18 1149 (1985). 155 M. Matsui, and M. Akaogi, Molecular Simulation 6 239 (1991). 156 A. Hallil, R. Tetot, F. Berthier, I. Braems, and J. Creuze, Physical Review B 73 165406 (2006). 157 R. Tetot, A. Hallil, J. Creuze, and I. Braems, Epl 83 (2008). 158 R. G. Hennig, T. J. Lenosky, D. R Trinkle, S. P. Rudin, and J. W. Wilkins, Physical Review B 78 (2008). 159 M. I. Baskes, and R. A. Johnson, Modelling and Simulation in Materials Science and Engineering 2 147 (1994). 160 G. Simmons, and H. Wang, Single Crystal Elastic Constants and Calc ulated Aggregate Properties (MIT, Cambridge, MA, 1971). 161 J. D. Gale, Journal of the Chemical Society Faraday Transactions 93 629 (1997). 162 J. D. Gale, and A. L. Rohl, Molecular Simulation 29 291 (2003). 163 J. A. S. Ikeda, Y. M. Chiang, A. J. Garrat treed, and J. B. Vandersande, J. Am. Ceram. Soc. 76 2447 (1993). 164 J. He, R. K. Behera, M. W. Finnis, X. Li, E. C. Dickey, S. R. Phillpot, and S. B. Sinnott, Acta Mater. 55 4325 (2007).
170 165 J. Nowotny, T. Bak, T. Burg, M. K. Nowotny, and L. R. Sheppard J. Phys. Chem. C 111 9769 (2007). 166 D. J. Wallis, N. D. Browning, P. D. Nellist, S. J. Pennycook, I. Majid, Y. Liu, and J. B. VanderSande, J. Am. Ceram. Soc. 80 499 (1997). 167 S. B. Sinnott, R. F. Wood, and S. J. Pennycook, Phys. Rev. B 61 15645 (2 000). 168 I. Dawson, P. D. Bristowe, M. H. Lee, M. C. Payne, M. D. Segall, and J. A. White, Phys. Rev. B 54 13727 (1996). 169 P. Y. Chiu, C. W. Lee, E. C. Dickey, S. R. Phillpot, and S. B. Sinnott, (submitted for review) (2011). 170 A. S. Barnard, and P. Zapol, Physical Review B 70 235403 (2004). 171 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). 172 M. Lazzeri, A. Vittadini, and A. Selloni, Physical Review B 63 155409 (2001). 173 M. Lazzeri, A. Vittadini, and A. Selloni, Physical Review B 65 119901 (2002). 174 M. Ramamoorthy, D. Vanderbilt, and R. D. Kingsmith, Physical Review B 49 16721 (1994). 175 P. J. Moller, and M. C. Wu, Surf. Sci. 224 265 (1989). 176 H. Onishi, and Y. Iwasawa, Surf. Sci. 313 L783 (1994). 177 P. W. Murray, N. G. Condon, and G. Thornton, Phys. Rev. B 51 10989 (1995). 178 R. J. Stokes, and D. F. Evans, Fundamentals of Interfacial Engineering (Wiley VCH, 1996). 179 U. Diebold, J. Leh man, T. Mahmoud, M. Kuhn, G. Leonardelli, W. Hebenstreit, M. Schmid, and P. Varga, Surf. Sci. 411 137 (1998). 180 J. Zhou, Y. C. Kang, and D. A. Chen, Surf. Sci. 537 L429 (2003). 181 U. Diebold, J. Lehman, T. Mahmoud, M. Kuhn, G. Leonardelli, W. Hebenstr eit, M. Schmid, and P. Varga, Surf. Sci. 441 137 (1998).
171 182 S. Fischer, A. W. Munz, K. D. Schierbaum, and W. Gpel, Surf. Sci. 337 17 (1995). 183 R. J. Stokes, and D. F. Evans, Fundamentals of Interfacial Engineering (Wiley VCH, 1997). 184 A. M. Asaduzz aman, and P. Kruger, Journal of Physical Chemistry C 112 19616 (2008). 185 W. J. Chun, Y. Koike, K. Ijima, K. Fujikawa, H. Ashima, M. Nomura, Y. Iwasawa, and K. Asakura, Chemical Physics Letters 433 345 (2007). 186 J. Graciani, J. J. Plata, J. F. Sanz, P Liu, and J. A. Rodriguez, Journal of Chemical Physics 132 (2010). 187 H. Aritani, N. Akasaka, T. Tanaka, T. Funabiki, S. Yoshida, H. Gotoh, and Y. Okamoto, Journal of the Chemical Society Faraday Transactions 92 2625 (1996). 188 H. Aritani, T. Tanaka, N Akasaka, T. Funabiki, S. Yoshida, H. Gotoh, and Y. Okamoto, Journal of Catalysis 168 412 (1997). 189 Y. Okamoto, H. Gotoh, K. Hishida, H. Aritani, T. Tanaka, and S. Yoshida, Applied Surface Science 121 509 (1997). 190 K. K. Bando, K. Sayama, H. Kusama K. Okabe, and H. Arakawa, Applied Catalysis, A: General 165 391 (1997). 191 N. Nomura, T. Tagawa, and S. Goto, Applied Catalysis, A: General 166 321 (1998). 192 D. Pillay, and G. S. Hwang, Journal of Molecular Structure Theochem 771 129 (2006). 193 U. Diebold, J. M. Pan, and T. E. Madey, Physical Review B 47 3868 (1993). 194 J. M. Pan, B. L. Maschhoff, U. Diebold, and T. E. Madey, Surface Science 291 381 (1993). 195 X. Z. Wu, R. Wang, and S. F. Wang, Applied Surface Science 256 3409 (2010). 196 W. R Tyson, and W. A. Miller, Surface Science 62 267 (1977).
172 197 D. G. Isaak, J. D. Carnes, O. L. Anderson, H. Cynn, and E. Hake, Physics and Chemistry of Minerals 26 31 (1998). 198 L. Koci, D. Y. Kim, J. S. de Almeida, M. Mattesini, E. Isaev, and R. Ahuja Journal of Physics Condensed Matter 20 (2008). 199 H. Z. Yao, L. Z. Ouyang, and W. Y. Ching, Journal of the American Ceramic Society 90 3194 (2007).
173 BIOGRAPHICAL SKETCH Tzu Ray Shan was born in Taipei, Taiwan, in 1978 to Wen Jing Shan and S h oo u Hw M aterials S cience and E ngineering from National Cheng Kung University, Tainan, Taiwan in 2001 and 2003 under the advisement of Profs. Kuan Zong Fang and Weng Sing Hwang respectively Following his commencement, Ray was qualified as a military officer in the Taiwan ( ROC ) Army and served mandatorily for 20 months as a platoon leader in Penghu Islands, Taiwan. After the military service, Ray worked as an assistant to patent agent/attorney for a year, d uring which time he fell in love with Ying Ju (Maggie) Chen who later became his wife in 2011 In August 2006, Ray was offered with the opportunity to pursue his Ph.D. degree at the University of Florida, and in May 2007, he was kindly provided with an as sistantship from Prof. Susan Sinnott and since then working under her advisement towards the degree. Ray earned his Ph.D. in M aterials S cience and E ngineering in May 2011.