|UFDC Home||myUFDC Home | Help|
This item has the following downloads:
1 DEVELOPMENT AND APPL ICATION OF REACTIVE POTENTIALS FOR THE ATOMISTIC SIMULATION OF MATERIAL SURFACES AND INTERFACES By BRYCE D. DEVINE 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 2010
2 2010 B ryce D. Devine
3 To Jocelyn
4 ACKNOWLEDGMENTS In an article describing potential development, P rof. Don Brenner states that this sort of work takes a person with peculiar perseverance. However, from my perspective, that is not necessarily true. What it takes is someone optimistically nave enough to begin in the first place and who is supported and nurtured by those with a peculiar sort of patience. I could not imagine complet ing this work and my studies without the constant and blessed support of my wife, Jocelyn. Through this time, I have had the encouragement of new solutions and to keep me motivated Joc elyn had no such encouragement. O nly a withering promise that I will finish this soon and yet her support has n ot wavered. In the end this work is ours together. Any achievement is Jocelyns as much as mine. Gratitude and accolades can not begin to recompense what she has given for me. Perhaps the best I can do is promise to never do any thing like this again. Along with Jocelyn, I also tried the patience of my advisor, Prof. Sinnott. This is not a trivial task for her patience is legendary. Her technical expertise, gift for teaching and constant encouragement were invaluable to this effort. I thank her sincerely for chance to work on these projects within this group. I will cherish and build upon the knowledge I have found here. I also had to good fortune to work with Prof Phillpot, whose expertise and experience were a constant source of new ideas and new directions Prof. S innott and Prof. Phillpot both set a lofty standard of scientific professionalism. The example they provide is inspirational and was as much an asset to this project as their guidance and insight. Several members of the University of Florida Computational Materials Science Focus Group have helped considerably in this work. Dr s. Alan McGaughey, Jia n guo Yu and Inkook Jang wrote many of the subroutines upon which I based my code. Coding is
5 an ungrateful job. They each spent years working through the intricacies of the codes It is labor for which the only acknowledgement they will receive is an esoteric listing in the source code. Ray Shan has been a tremendous part ner in developing these models. With out our daily discussions, I doubt these projects would be finished. Some of this work also could not be finished without the assistance and expertise of D onghw a Lee and Travis Kemper and the suggestions of Dr. Tao Liang to all of whom I am grateful. Lastly, funding for much of this work was graciously provided by the National Science Foundation and the Department of Defense. All work was performed on the computers operated by the Computational Materials Science Focus Gr oup at the University of Florida. With out this resource and the dedicated graduated student who maintain it, we would accomplish very little.
6 TABLE OF CONTENTS Page ACKNOWLEDGMENTS .................................................................................................. 4 LIST OF TABLES ............................................................................................................ 8 LIST OF FIGURES ........................................................................................................ 10 ABSTRACT ................................................................................................................... 11 1 ANALTYICAL POTENTIAL S IN COMPUTATIONAL M ATERIAL SCIENCE ........... 14 1.1 The Promise of Multiphase Atomistic Simulations ............................................ 14 1.2 Analytical Potentials for Multi Material and Multi phase Simulations ................ 16 1.3 Variable Charge Reactive Potentials ................................................................ 18 1.4 An Overview of the Dissertation ........................................................................ 20 2 THE THEORETICAL FOUNDATION OF CHARGE OPT IMIZED MANY BODY POTENTIALS ......................................................................................................... 25 2.1 The Quantum Mechani cs of Inter Atomic Interactions ...................................... 25 2.2 Wavefunction Theory (WFT) ............................................................................. 29 2.3 Density Functional Theory (DFT) of Inter Atomic Interact ions .......................... 48 2.4 A Self Consistent Classical Electrostatic Potential ............................................ 57 3 A MULTILEVEL COMPUTATIONAL ANALYSIS OF FLUOROCARBON POLYATO MIC DEPOSITION ON DIAMOND ......................................................... 85 3.1 The Computational Study of Polyatomic Deposition Reactions ........................ 85 3.2 Computational Details ....................................................................................... 90 3.3 Simulation Results ............................................................................................ 95 3.4 Discussion of Results ...................................................................................... 103 3.5 Concludi ng Remarks ....................................................................................... 109 4 METHODOLOGY FOR ATOM ISTIC SIMULATIONS OF COPPER OXIDATION USING A CHARGE OPTIM IZED MANYBODY (COMB) POTENTIAL ................. 116 4.1 Current Theoretical Models of Copper Oxidation ............................................ 116 4.2 Computational Methodology ........................................................................... 119 4.3 Potential Parameterization .............................................................................. 121 4.4 Results and Discussion ................................................................................... 125 4.5 Concluding Remarks ....................................................................................... 136 5 CHARG E OPTIMIZED POTENTIALS FOR ALUMINUM AND ALUMINUM OXIDE SYSTEMS: SUCCESSES, CHALLENGES AND OPPORTUNITIES ........ 149
7 5.1 Background and Motivation ............................................................................ 149 5.2 Challenges with Current Theoretical Models for Alumina ............................... 150 5.2 Parameterization of the COMB Potential ........................................................ 155 5.5 Results and Discussion ................................................................................... 159 5.5 Concluding Remarks ....................................................................................... 165 6 CURRENT CAPABILITIES AND FUTURE POSSIBILI TIES OF CHARGE OPTIMZED MANY BO DY POTENTIALS .............................................................. 179 6.1 The COMB Code ............................................................................................ 183 6.2 The Fitting Problem ......................................................................................... 185 6.3 Future Applications ......................................................................................... 186 LIST OF REFERENCES ............................................................................................. 188 BIOGRAPHICAL SKETCH .......................................................................................... 198
8 LIST OF TABLES Table page 3 1 Geometric parameters for CF23 radicals and cations. ...................................... 113 3 2 Classical MD predicted reactions between C F3 and diamond (111) ................. 113 3 3 DFT MD predicted reactions between CF3 and the diamond (111) ................. 113 3 4 Reaction enthalpies (in eV) between CF3 and adamantane (C10) ..................... 114 3 5 Reaction enthalpies (in eV) between CF3 and C22 ............................................ 115 3 6 R eaction enthalpies (in eV) bet ween CF3 and a periodic diamond slab. .......... 115 4 1 Atomic and electrostatic potential parameters .................................................. 144 4 2 Bond dependent potential parameters .............................................................. 144 4 3 Coefficients for Legendre polynomials .............................................................. 145 4 4 Properties of metallic copper ............................................................................ 145 4 5 Ionization potentials and electron affinities (in eV) for Cu and O ...................... 146 4 6 Properties of m olecular o xygen ........................................................................ 146 4 7 Properties of Cu2O (Cuprous Oxide) ................................................................ 147 4 8 Properties of m onoclinic CuO ........................................................................... 147 4 9 Hf (in eV) of c opper oxides ............................................................................. 148 4 10 HAds (in eV) of O2 on Cu(100) surface ............................................................ 148 5 1 Atomic and electrostatic potential parameters .................................................. 174 5 2 Interaction dependent potential parameters ..................................................... 174 5 3 Coefficients for Legendre polynomials .............................................................. 175 5 4 Properties of metallic aluminum ........................................................................ 176 5 5 Ionization potentials and electron affinities (in eV) for Al and O ....................... 177 5 6 Properties of Al2O3 ........................................................................................ 177
9 5 7 Interfacial separation distance, Do (in ) and Wsep (in J.m2). ............................ 178
10 LIST OF FIGURES Figure page 1 1 Si nanocluster in amorphous SiO2. ................................................................. 23 1 2 The Modelers Toolbox ..................................................................................... 24 2 1 The interaction energy between two charged atoms. ......................................... 83 2 2 Self energy functions for O and Cu. .................................................................... 84 3 1 Target sy stems are used for calculation of reaction enthalpies ........................ 111 3 2 Impact sites on the diamond(111) surface. ....................................................... 112 3 3 The geometry of the transition state ................................................................. 112 4 1 Crystal structures for the two low energy oxide phases .................................... 137 4 2 Surface charge in metallic copper.. .................................................................. 137 4 3 The bond dissociation energy of O2 anions ..................................................... 138 4 4 The energy per unit volume of the three low energy phases ............................ 139 4 5 O2 is placed on the Cu(100) surface in the six orientations .............................. 140 4 6 The relaxed reconstructed Cu(100) surface. .................................................... 141 4 7 The Cu2O(111)||Cu(100) interface. ................................................................... 142 4 8 Planar average charge density ........................................................................ 143 5 1 The crystal structure of Al2O3 (A) .................................................................. 167 5 2 The Al O Al bond angle distribution in corundum ............................................. 168 5 3 The energy vs. volume of several aluminum oxide phas es. ........................... 169 5 4 A closed packed plane of Al atoms ................................................................... 170 5 5 A second layer of Al atoms in the FCC metal ................................................... 171 5 6 The interfaces with the highest (A) and lowest (B) values ............................... 172 5 7 The charge density (blue line) and separation distance ................................... 173
11 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 APPL ICATION OF REACTIVE POTENTIALS FOR THE ATOMISTIC SIMULATION OF MATERIAL SURFACES AND INTERFACES By Bryce D. Devine December 2010 Chair: Susan Sinnott Major: Materials Science and Engineering Interest in atomic scale computational simulations of multiphase systems has grown substantially i n recent years as our ability to simulate multi million atom systems has become commonplace. The main limitation is now the lack of a theoretical framework that can smoothly model atoms in dissimilar bonding conditions such as across a metal metal oxide interface. Recently developed charge optimized many body potential (COMB) have addressed several of the technical challenges to atomistic simulations of systems composed of multiple phase and multiple materials interacting together. The COMB formalism merges variable charge electrostatic interactions with a bond order potential that has the capacity to adaptively model metallic, covalent and ionic bonding in the same simulation cell. Presented here is the development of an inter atomic COMB potential ener gy function from its basis in quantum mechanics up through its adaptation and parameterization for two challenging simulation problems: the oxidation and oxide island formation on copper surfaces and the aluminum aluminum oxide interface. In theory, COMB is an analytical extension of tight binding theory to a fully classical determination of inter atomic interactions. In practice, it is
12 found that several empirically derived functions are also necessary to adapt the potential to systems of interest. In this work, these additional functions are chosen to replicate interactions that are lost in the approximations made in the derivation of the COMB potential. The resulting COMB formalism accurately reproduces key observables in both materials such as mechanic al properties, phase order, formation enthalpies, surface energies, and defect formation enthalpies. The performance of the potentials indicates that the COMB potential holds significant promise in extending our abilities to simulate multiphase systems. M olecular dynamics simulations of the oxidation of the Cu (100) surface using the COMB potential reveal a fairly low occurrence of reactions. Single point calculations of O2 on this surface predict that molecular dissociation only occurs when both O atoms are close to adjacent four fold hollow sites. The molecule dissociates to form a c(2x2) surface at low temperature and low oxygen coverage. At room temperature and with higher coverage, the surface reconstructs to a missing row 45 ) 2 2 2 ( R co nfiguration. These findings are consistent with experimental observations. Similarly COMB based molecular dynamics simulations of variations of crystalline Al(111)|Al2O3(0001) interfaces predict atomically sharp interfaces that are stable through the melt ing of the Al metal. The potential does not clearly distinguish differences in interfacial adhesion with different orientations of the metal to the oxide. This is inconsistent with first principles calculations and warrants further investigation. One consideration with advances in computational power and improved algorithms is first principles based methods may be adequate to model multiphase systems. This possibility is explored through a multilevel study of fluorocarbon
13 deposition on diamond using bot h classical molecular dynamics (MD) using the reactive empirical bond order potential (REBO) and density functional theory based MD (DFTMD). The comparison highlights that classical methods are still several orders of magnitude more efficient than current DFT methods, thus yielding better statistics, but DFT models are more reliable in their predictions. In this study, it is found that CF3 and CF3 + react differently when deposited with energy of 50 eV. It is also found that CF3 can add directly to the hy drogen terminated diamond surface when coupled with the simultaneous evolution of HF. The single step addition reaction is only predicted with DFT due to short interaction distances in REBO. The results strongly suggest an improved REBO potential that inc ludes charge optimized electrostatic interactions and better represents the long interaction distances in key transition states would greatly benefit our understanding of this commercially relevant process.
14 CHAPTER 1 ANALTYICAL POTENTIAL S IN COMPUTATIONAL MATERIAL SCIENCE 1.1 The Promise of Multiphase Atomistic Simulations As engineering progresses at the nanometer length scale, computational simulations are finding an increasingly significant role in the design and optimization of new electronic and mechanical devices. This is due in part to the confluence of two favorable trends. First, as the size of devices decreases, the total number of atoms in the structure is reduced. In a nanoscale device, the total number of atoms is on the order of tens of mi llions to billions. At the same time, with the relentlessly consistent year over year increase in computer performance, it now fairly routine to conduct entirely atomistic simulations on systems composed of tens on millions of atoms. Multibillion atom sim ulations are possible on the current generation of the most powerful supercomputers. If we hold to our current course of progress, in a relatively short time life size simulations of entire nanodevices will become routine rather than leading edge. While increases in computer power expands the accessible length scale across all levels of theory, the bulk of computational work for lifesize systems defaults to empirical atomistic methods in which interactions between all atoms in the systems are explicitly described with classical Newtonian mechanics. Such methods are highly efficient yet provide atomic scale detail regarding the structure and energies of critical features, such as interfaces, and the atomic scale processes that affect material performance t hroughout the device. In reality, inter atomic interactions are best described with quantum mechanics based methods that consider all electrons in the systems. Quantum based methods are generally regarded to be more accurate with a
15 higher degree of reliab ility in their predictions than atomistic potential. However, electronic structure calculations are limited by their computational demand to systems of just a few thousand atoms. At the other end of the length scale in materials simulations are mesoscale methods which look at materials properties at the microscale. Such methods lack the atomic scale detail that is highly valuable to materials selection and optimization at the nanoscale. To illustrate this point, Figure 1 1 shows an image from an atomi stic simulation of a silicon nanocluster embedded in an amorphous silica matrix. The colors correspond to the coordination of the individual atoms. Normally coordinated atoms in the silica matrix are not visible for clarity. It has been discovered recently that systems with nanoclusters such as this exhibit several intriguing optical properties including the ability to emit light13. While the commercial significance of a source of solid state lighting made from some of the cheapest and must abundant materials on earth is easily appreciated, the technological barriers to developing a commercially viable device are huge and quite daunting. As it is currently, there is no clear agreement on how emission occurs in the material. One school of thought proposes quantum confinement effects within the cluster or at the surface give rise to emission46, while others attribute it to under coordinated defects at the interface and within the amorphous matrix79. It is nearly impossible to answer the question solely through experimentation. There is currently no reliable experimental technique that can resolve phenomena at buried interfaces at the atomic scale in situ. Unfortunately, the system is quite challenging for theoretical methods as well. The nanocluster in Figure 1 1 is 4 nm in diameter which is a realistic size for the material. The simulation cell contains ~60,000 atoms between the cluster
16 and the matrix. The system size is minimal for a cluster of that size, yet the system is nearly 6 times larger then what can be studied through quantum based methods. Continuum based models lack the atomic scale detail necessary to answer the critical questions such as what is the structure and coordination at the interface. The ideal method to study this system is an atomistic analytical potential designed to model both phases simultaneously. Naturally, life size atomistic simulations of any device most likely involve multiple interfaces. In most cases, such simulations include a surface, which is at least an interface between a solid and a vacuum. More realistically, devices are multifunctional and involve the integration of several materials that contribute different functionality to the overall operation. Therein lays the crux with the present state of theart of atomistic simulations. Currently, what is lacking is a theoretical methodology that can simultaneously model multiple phases and all types of chemical bonding: ionic, covalent, metallic and van der Waals interactions. 1.2 Analytical Potentials for Multi Material and Multi phase Simulations In this dissertation, the classical methodology that describes inter atomic interactions is called a potential as is the convention of materials scientist s and solid state physicists; in molecular systems, the methodology is conventionally referred to as a force field. While describing molecular interactions is a goal of this work, the majority of the systems we seek to model are in the solid state and so I will stick with the conventions of the solidstate community. Typically a potential is designed and parameterized to describe one type of chemical bonding. Over time, this has lead to several categorizations of potentials based on the type of material for which they are applicable: Buckingham for ionic
17 materials ,10 Tersoff potentials for covalent bonding11, 12 or the Embedded Atom Method (EAM) for metallic systems13 are well documented examples. The obvious limitation with this scheme of development is that none of these methods can seamlessly model several bonding environments, such as what occurs in a heterogeneous interface between dissimilar materials. Recently, this limitation has been breached by several advances in potential development. While the goal of a universal method applicable to all bond types has not yet been realized, significant progress has been made based on two key developments. The first is a rather new class of models called variable charge potentials in which the partial charges of individual atoms are not fixed but instead are determined in a self consistent m anner based on the principle of electronegativity equilibration (Qeq).14 Essentially, variable charge potentials describe the electrostatic interactions between atoms as a balance between the Coulombic interaction between charge and the cost to fro m a charge on an atom. This allows for the simultaneous handling of an atom in various oxidation states within the same classical simulation. Examples in the literature can be found that apply Qeq type electrostatic schemes to all bond types depending on the type of short range potential with which the scheme is coupled.1520 The Qeq methods describe only the electrostatic portion of the total energy in a system and must be added to a nonelectrostatic potential to describe short range contributions to the total energy. For example, Streitz and Mintmire coupled a Qeq scheme to a Finnis Sinclair21 potential for metals to model aluminum/alumina interfaces.22 Recently, bond order type potentials23 and especially reactive bond order potentials16, 24, 25 that describe bond breakin g and new bond formation in covalent
18 materials, have been shown to have the flexibility to model several bond types from metallic to purely covalent with varying degrees of covalency in between. This is demonstrated by the work of Iwaski15 and later Yu et al.18 who applied a Tersoff type bond order potential to various transition metals. Coupling of a Qeq electrostatic scheme with a bond order potential provides a means to extend the method to ionic type bonding. Such an approach was developed by Yasukawa17 who paired Qeq with a Tersoff potential. More recently Yu et al. developed their charge optimized many body (COMB) formalism based on the work of Yasukawa.18 The initial version of COMB fixed a few instabilities in the Yasukawas potential. More recently, Shan et al. corrected a few of the i nstabilities in the original COMB formalism.26 In particular, they parameterized a version of the potential for hafnium and the various phases of hafnium oxide, which, due to their complex crystallography, are difficult materials to simulate.26 Recently an article by Phillpot and Sinnott discussed the promise of variable charge bond order potentials for multiphase simulations.26 The capabilities that this approach offers materials modelers is illustrated in Figure 1 2, which is published in that article. The figure presents the three main bonding types to which most current potentials are limited: metallic, ionic and covalent, and the technological problems we can address with a multiphase atomistic potential. 1 .3 Variable Charge Reactive Potentials In addition to the extensibility to new materials systems, variable charge potentials offer several advantages over traditional atomistic potentials in cases where bonds are broken and formed. Currently the main appli cation of reactive potential is for hydrocarbon system where the reactive empirical bond order potential (REBO) has been used extensively.24, 25 REBO is an extension of the Tersoff potential11, 12 designed
19 to model bond breaking in hydrocarbon systems. The main limitation of REBO is its lack of a description for electrostatic interactions which restricts its use to charge neutral covalently bound systems. The Reax FF force field is another reactive potential based on an extended Tersoff formalism that was originally designed for hydrocarbon systems.16 ReaxFF includes variable charge electrostatics and has been extended to both metallic and ionic systems.19 The main limitation with ReaxFF is its scalability, which limits its practical application to systems with just a few thousand atoms. However, the success of ReaxFF in modeling mixed bond states further supports the concept of combing variable charge electrostatics with a bond order potential to develop a multiphase potential. Reactive potentials are often noted for their efficiency which is often the first consideration in choosing an atomistic level of theory. Variable charge techniques add several key characteristics to extend the applicability of reactive potentials, the first and foremost of which is the ability to model various oxidation states. Take, for example, an ionic system described with a fixed charge model. In such a scheme, each atom has a set charge regardless of the local environment in which it exists. The charge is the same whether an atom is bound in the solid, at the surface or liberated to infinite distance. In a fixed charged ionic model, the energy to completely atomize the system is actually the ionization energy defined as the energy required for break the material into infinitely separated ions. The more chemically meaningful measure is the enthalpy of formation, which is the energy difference between the final material and its constitutive elements in their standard reference states. The enthalpy of formation is unobtainable from a fixed charge potential since the standard reference states for
20 elements are charge neutral. A variable charge potential is the only atomistic method capable of calculating formation enthalpies relative to the standard reference state in ionic materials. The other main application for empirical methods is in cases where higher fidelity methods are not applicable. In computational materials science, the most frequently applied quantum mechanics based methodology is KohnSham density functional theory (DFT) ,27, 28 which is especially well suited for metals and heavier atoms i n the solid state. In several instances, we seek to look at a problem for which DFT is prone to error. Bond dissociation is a typical example where static correlation effects, which DFT misses 29, contribute significantly to the total energy at extended bond lengths. An empirical potential can be parameterized to a much higher level of theory that captures such effects and efficiently model larger systems that are not tractable by any other means. 1.4 An Overview of the Dissertation The most reliable analytical potentials are derived from the quantum physics that govern chemical bonding. This is the case for most of the generally used potentials for materials simulations such as EAM and Tersoffs bond order. While it is possible to produce a potential based on empirical observations, chance are good that it will have limited applicability beyond what it was parameterized for. The second chapter reviews the quantum physics that describe chemical bonding as it pertains to potential development. A COMB potential formalism is derived through classical approximations to the quantum physics based interatomic interactions. The intent to illustrate the connectivity between the derived classical potential energy function and DFT.
21 Chapter 3 of this work compares the application of a reactive empirical potential and DFT to the simulation of the deposition reactions on diamond. This is an area of study where atomistic simulations have provided significant ins ight into the reactions the govern deposition products. The work highlights the strengths and limitations of the various methods and supports the further development of COMB type potentials. The Chapters 4 and 5 cover the development and application of COMB potentials for two relevant materials problems. Chapter 4 presents a potential designed to model the oxidation of copper and the formation and growth of copper oxide islands on the metallic surface. Metal oxidation is a commercially and technological ly significant process for which copper has been a model system .3034 More recently, it has been found that the size and shape of nanometer sized oxide islands grown on the copper surface can be influenced by the conditions under which they form .3537 The overall goal is that with a bit more insight into the mechanism and physics that govern island formation we can find a means to dictate the size in and shape of these nano scaled structures. To date, the limitation for simulations of these systems has been the lack of an atomistic method to apply to the problem. Chapter 5 presents a potential for the simulation of aluminum/alumina interfaces. Due to its technological significance, this material system has been the subject of several development effor ts.36, 38, 39 Many potentials fail t o predict the correct ground state crystal structure of alumina, which limits their reliability in predicting interfacial structures. There are currently two potentials that predict the corundum structure as the ground state. The ReaxFF formalism has been adapted to the aluminum and alumina.40 The s calability limitations of ReaxFF limit its use to small interfacial systems. The other
22 model is the compressible ion model (CIM)41, which includes the polarization effects that are believed to stabilize the corundum structure.39 However, the model does not include variable charges, whic h precludes its use from interfacial simulations. The Chapter 6 summarizes the current state of the art COMB potentials and the areas for future development. This is an ongoing and active are research with several projects seeking to extend the parameteri zation to more complicated system, such as hydrocarbons and rare earth oxides. The intent is to close with a summary of finding and a review of the tools that have been developed to assist in these efforts.
23 Figure 1 1. Si n ano cluster in amorphous SiO 2 This snapshot from an atomistic simulation shows the interface Si nanocluster embedded in amorphous SiO2. Atoms are colored based on coordination number (CN). Normally coordinated atoms are not visible. CN=2 CN=3 CN=1
24 Figure 1 2. The Modelers Toolbox as describe by Sinnott and Phillpot in a recently published article discussing simulations of multiple material systems. The image illustrates the technical problems that can be addressed with a potential, such as COMB, that has the capability of modeling more than one bond type. Reprinted with permission from S. R. Phillpot and S. B. Sinnott, Science 325 1634 (2009).
25 CHAPTER 2 THE THEORETICAL FOUNDATION OF CHARGE OPT IMIZED MANY BODY POTENTIALS 2.1 The Quantum Mechanics of Inter A tomic Interactions A main goal of the work with COMB is to develop a potential formalism with the flexibility to model dissimilar materials and with the transferability to describe the varying bonding environments simultaneously. The most promising approach to accomplish that objective is to begin with a potential that is derived by making classical approximations to the quantum physics involved in inter atomic interactions This is essentially how Abell and Pettifor derived their respective analytical bond order potentials from tight binding density functional theory (TB DFT).23, 42 This style of potential is perhaps the most flexible formalism in common use today and has been adapted with slight modification to covalently bound semiconductors ,11, 36 hydrocarbons ,16, 24, 25, 43 ionic material19, 44 and metals .15, 18 Tersoffs version of the Abell formalism is one of the most reliable potentials for Si and other semiconductors.11, 12 Brenner built upon the Abell Tersoff formalism to develop the Reactive Bond Order Potential (REBO)24 and its second generation25 for hydrocarbon systems. The reactive term in the name indicates that this potential has the flexibility to model bond breaking in carbon systems. Similarly, the ReaxFF force field, which too is based on Tersoffs model, has been extended to a wide range of materials spanning nearly all regions of the periodic table.16 Pettifors analytical bond order formalism was originally developed to distinguish single double and triple bonds as well as radicals in carbon materials42 but has since been adapted to more complex covalently bound systems such as GaAs and SiC and metals .35, 36, 4548 Brenner has shown how the embedded atom methods (EAM) originally designed for metallic systems are algebraically equivalent to Abells model.13
26 Recent versions of EAM extend this formalism further to covalent systems such as Si and C .49, 50 All together, the analytical bond order formalism and its derivatives have been adapted to the most comprehensive list of materials of any current potential formalism. While the analytical bond order form alism has demonstrated the flexibility to model many systems, historically it has lacked the transferability to simultaneously model several phases. Variable charge schemes address this limitation. The flexibility of the bond order formalism is attributabl e to its theoretical origins in DFT.21, 51 A dding additional electrostatic interactions poses the risk that the additional interactions may compromise the flexibility of the formalism. However, as with the bond order formalism, variable charge schemes are shown here to have a theoretical basis in DFT as well. The coupling of the two methods can be viewed as incorporating a few of the contributions to the bond energy t hat are disregarded in bond order formalism. In a review of the development of the REBO potential, Brenner proposes that transferability and flexibility are best achieved when a potential replicates the essential quantum physics of chemical bonding.52 Finnis emphas izes this point in a similar review of Pettifors bond order formalism and in his book where he derives many of the most effective interatomic potentials from TB DFT.21 The coupling of a variable charge electrostatic potential with a bond order potential follows this same approach. When developing a potential, it is always tempting to add additional functions i n an ad hoc manner to correct aberrant predictions or improve the behavior of the potential. While corrections are in most cases unavoidable, they often have unintended consequences and must be chosen with care. In his review, Brenner suggests four
27 charac teristics of effective potentials that must be balanced while developing the functional form of the potential: Flexibility: The potential should be able to reproduce a wide range of materials properties such as cohesive energies, surface energies, defect formation energies, etc. Accuracy : The potential should reproduce the properties in the fitting database within the designed tolerance. Transferability : The parameter set should be applicable to a wide range of bonding environments and other phases of the material. Computational Efficiency: The appeal of empirical potentials over more transferable first principles methods is the speed of the calculations and scalability to larger systems. There is a lower bound on the efficiency of an empirical potential where, if it is nearly as costly as first principles calculation, there is little reason to spend the effort fitting and developing it. When we add functions to improve performance of the potential for one phase we naturally risk compromising the transferability to other phases. This especially pertinent when designing potentials for interfaces where the parameters must be transferable to several phases in a broad range of coordination environments. Take for example, a graded interface between an oxide and a closepacked metal. In the interface region, coordination ranges from 12 nearest neighbors in the pure metal to substantially less in the stoichiometric oxide. Between these boundaries, in the interface region, there is no constraint on stoichiom etry. Furthermore, variable charge electrostatics also lifts any constraint on charge neutrality at the interface; the system as a whole can be constrained to be charge neutral, but in the interface region, nonneutrally charged local
28 phases may develop. Achieving transferability to all these unusual and unexpected possible phases is a delicate task that requires well chosen functional forms. Success is most likely when the additional functions are chosen to replicate the interactions as described by the underlying quantum physics. Transferability of a potential also depends on the fitting process. A tempting assumption is that transferability can be achieved solely by fitting to a comprehensive database. While fitting is a critical part of the potential development process, where accuracy and the flexibility to describe a wide range of phases are necessary for reliable results, it is often quite possible to achieve an accurate fit to properties that should be indistinguishable by the potential. A typical example is the ability of an EAM potential to stabilize a face centered cubic (FCC) crystal structure relative to a hexagonal close packed (HCP) arrangement. The physics that an EAM potential describes cannot distinguish between these phases. However, w ith careful fitting of the cutoff distance, the FCC structure can be stabilized relative to HCP. In such a case, transferability is achieved by fitting a convenient cutoff that selects the interactions that give a decent result rather than accurately desc ribing the energy per interaction. It is doubtful that such a potential will be transferable to the unusual bonding situations that arise at surfaces and interfaces or to the metal metal interactions in an oxide. The point here is that while fitting is a critical phase of the development, it is possible to fit the correct behavior for the wrong reasons. Transferability is best achieved by beginning with a well designed functional form that replicates the underlying physics and fitting that potential to a w ell chosen data set.
29 2.2 Wavefunction Theory (WFT) 2.2.1 The time independent Schr dinger equation The simplest description of a material at the quantum level is as a system of nonrelativistic electrons acting in response to an external potential Vext com prised of the positively charged nuclei. One way to describe the physics of the electrons is through wave mechanics, where properties of the chemical system ar e described by a wavefunction.53 Observable properties of the system, such as energy, are found when an appropriately chosen operator acts upon the wavefunction. In the case of nonrelativistic electrons, the energy of the system is given by the equation: E H (2 1) This is the time independent Schr dinger equation in its most basic form. Here E is the total energy of the system and the wavefunction, ( ri, si) is an eigenfunction of position and spin for each electron in the system. Equation 21 is an eigenvalue equation where an operator, in this case, acts on the wavefunction to yield a scalar eigenvalue, E The Hamiltonian operator in Equation 21, is the operator that gives the total energy for the system. When written out explicitly the Hamiltonian is: I i J J I J I i i j j i i I i i I I I I i i eR R Z Z e r r e R r Z e m m H0 2 0 2 0 2 2 2 2 24 2 1 4 2 1 4 2 2 (2 2) where, i and j sum over electrons located at position ri (xi,yi,zi) I and J sum over nuclei at positions RI, me is the mass of an electron, mI is the mass of the nuclei, Z is the atomic number, e is the charge of one electron, 0 is the vacuum permittivity. The term
30 is Plancks constant divided by 2 and the term 2 i s the Laplacian operator defined in Cartesian coordinates as: 2 2 2 2 2 2 2 i i i iz y x (2 3) The energy operator contains components for both the kinetic and potential energy of the electron and nuclei. The last three terms in Equation 22 describe the potential energy and are essentially the same as the classical electrostatic interactions. The factor of in the last two terms accounts for double counting since the summations run over each interaction twice. At the quantum level, the kinetic energy is not equivalent to its classical analog, but, instead is determined as the eigenvalue of the kinetic energy operator54: 2 2 2 m T (2 4) The wave function, is an obtuse concept. Essentially, it is a mathematical function whose arguments are the positions and spins (r,s ) of the electrons. The s pin, si, introduces two components for the description of each electron, conventionally referred to as (si) and (si) or rather the spin up and spin down component for each electron. The main effect from spin is that the wavefunction is anti symmetri c, meaning the wavefunction changes sign when the positions of two electrons are exchanged.55, 56 The product of the wavefunction with its complex conjugate, (*) represents a probability density. As an eigenfunction, could be multiplied by any scalar quantity and still return a true solution, thus, giving an infinite number of solutions to Equation 21. In practice, the wavefunction is crafted to be orthonormal, meaning the components
31 are orthogonal to one another and the wavefunction is multiplied by a unique normalization constant that satisfies the condition: ij j idr (2 5) w here, ij is the Kronecker delta function, which is equal to one w hen i=j and zero otherwise. Equation 25 uses shorthand notation where dr indicates integration over all spatial components and summation of the spin components, s. The normalization of the wavefunction leads to a representation of the electron density at any position, ri, as: dr N ri 2) ( (2 6) N in Equation 26 is the number of electrons. 2.2.2 The BornOppenheimer approximation The complexity Equation 22 can be reduced with a few key approximations. First, since in most circumstances the mass of the nuclei are significantly larger than the mass of an electron and the electrons move at a significantly higher velocity than the nuclei, the motion of the nuclei may be neglected when determining the energy of the electrons. This approximation, which is referred to as the BornOppenheimer approximation, drops the second and fifth terms from Equation 22. With this approximation, the equation yields a total electronic energy, Eel, for the system rather than the actual total energy for the system .57 A fundamental hypothesis of quantum mechanics is that there has to exist a lowest energy state for the system called the ground state, Eo, which is described with the corresponding ground state wavefunction, o. 58 If this were not true then an infinitely low energy could be achieved by simply packing the electrons in closer to the nuclei while maintaining the kinetic energy at a constant value. A related approximation, the
32 adiabatic approximation, assumes that under normal conditions, on the time scale of nuclear motion, electrons relax to the ground state instantaneously.59 Under the adiabatic approximation, the energy of any arrangement of nuclei can be determined with the ground state configuration of electrons. The two approximations introduce inconsequential error in most circumstances, yet have several significant ramifications in regards to empirical potentials. Together the two approximations allows us to clearly separate electronic and nuclear motions, which, in turn, allows us to derive an empirical potential for inter atomic forces based upon the ground state configuration of the electronic structure for any arrangement of nuclei. For any given set of nuclear positions, there is one unique ground state electronic configuration. Consequently a finite difference in energy can be determined for various stresses and perturbations on the system. Without the approximations, any derived potential would be based on a one or several most probable configurations rather than the unique ground state. The appropriateness of approximation is contingent of course on the definition of normal conditions. In most applications we are concerned with low energy processes around standard conditions. The BornOppenheimer approximation sets an upper bound in energies below which it is valid. Processes that involve high nuclear velocities or sufficiently high energies to allow excited electronic states to participate are somewhat spurious. This limitation carries through to any empirical potential derived from the ground state configuration. For most applications involving static calculations or low energy dynamics around standard temperature, the approximation is very reasonable. In computational materials we can run into difficulties when we look at very
33 high energy processes such as high energy depositions or radiation bombardment in which case results may be unreliable. The separation of nuclear and electronic motions allows us to simply the expression of the Hamiltonian in Equation 22. First we can drop all nuclear nuclear interactions from the calculation. Secondly, we can define the potential generated by the field of the positively charged nuclei and all other potentials t hat are not dependent upon the electrons as an external potential, Vext: I I I extR r Z r V ) ( (2 7) The external potential and any other potential acting on the electrons can be defined as an operator that gives the expectation value of the electronn uclear energy contribution to the total electronic energy of the system. The expectation value for the electronnuclear Coulombic attraction is determined by operating extV on the wavefunction, premultiplying by the complex conjugate of the wavefunction and integrating over all special components: dr V Eext eZ (2 8) The interaction can also be described as the interaction between the Vext(r) and electron density as defined in Equation 26 which will be used later: dr r V r Eext eZ ) ( ) ( (2 9) 2.2.3 The variational principle If the correct wavefunction is known, then any observable value of the system can be determined by operating on the wavefunction with the appropriate operator. The trick is finding the appropriate wavefunction for the system. To so do, we can exploit
34 the idea that there exists a lowest energy ground state, which sets a lower bound on any trial solutions. The best solution, therefore, corresponds to the lowest energy, which set up a constrained minimiz ation problem to determine the ground state configuration. This concept, called the variational principle,60 is an essential concept in most quantum mechanical calculations and is well described in the standard text books. I will use the proof by Cramer to illustrate.58 A normalized trial wavefunction, can be expressed as a sum of basis functions multiplied by an expansion coefficient, Ci. As an example, a trial wavefunction can be constructed as a linear combination of orthonormal single particle wavefunctions, I: i i iC. (2 10) Since is normalized, the basis functions and coefficients are constrained to satisfy the orthonormal condition meaning i i ij ji j i j i ji j ic c c dr c c dr2 21 (2 11) To find the energy of any particular state, both sides of Equation 22 can be multiplied with a trial wavefunction and integrated over all spatial components: dr E dr Hi (2 12) Since E is a scalar, it can be moved out of the integral. Using the orthonormal condition, Equation 212 can be recast in terms of the coefficients, Ci
35 i i i ij i ij j i i j i ij j i i j ij j iE c E c c dr E c c dr H c c dr H 2 (2 13) The results of Equations 211 and 213 can be used to construct an expression for the energy of any state relative to the ground state as: ij i iE E c dr E dr H ) ( 0 2 2 0. (2 14) The left hand side of Equation 214 is always greater than or equal to zero since c2 0 for all real values of c, and (Ei E0) 0 b y the definition of E0. This sets up the following inequality: 0 2 0 dr E dr H. (2 15) With rearrangement, Equation 215 becomes : 0 2 E dr dr H (2 16) In the example above, the wavefunction describes a single particle, but this is not a limitation of the method. The variational principle applies to single particle wave functions as well as many electron wavefunctions describing many atom systems. Also, in this example, the trial wavefunction was composed of a linear combination of orthonormal wavefunctions. This too is not a constraint of the principle but rather a convenient condition chosen for illustrative purposes. 2.2.4 Antisymmetry and determinant wavefunctions In the example provided in Section 2.2.3, the trial wavefunction is composed of a linear combination of basis functions that approximates the actual wave function. A
36 similar approach can be used to construct a many electron wavefunction. The choice of a trial wavefunction must consider that each electron has both spatial and spin components. A trial a wavefunction, (x1,xn) that captures both components can be constructed as linear combination of single particle wavefunctions a(x1) composed as a product of spatial and spin components: ) ( ) ( ) (i i a i as r x (2 17) In Equation 2 17, the arguments xi refer to spinorbitals that include both spatial ,a(r) and spin components, (s) .The spin component arises as a consequence of Diracs relativistic quantum mechanics61 and result in the wavefunction being antisymmetric, meaning the wavefunction changes sign with the interchange of any electron pair: ) ( ) (j k i k j ix x x x x x (2 18) The antisymmetric nature of the wavefunction also leads to Paulis exclusion principle, which states, in basic terms, that no two electrons can be described with the same identical quantum numbers.62 A first choice to consider for a trial wavefunction would be product of spin and spatial components. To preserve orthonormality in ( x) the spin component s must also satisfy the conditions: 0 ) ( ) ( 1 ) ( ) ( ) ( ) ( i i i i i is s s s s s (2 19) Normalization must also consider integrating over all spin orbital states: 11 2 ndx dx (2 20) The electron density when determined according to Equation 220 must also include the spin components such that the density at position r1 is :
37 1 ) (2 1 2 1 ndx dx ds N r (2 21) An effective means of maintaining antisymmetry in the wavefunction was developed by Slater who represented the many electron wavefunction as determinant.63 T his approach exploits the property of determinants that the sign of the determinant changes when two rows or two columns are exchanged. For example the wavefunction for a system composed of two electrons with the same spin in different single particle sta tes can be represented as: ) ( ) ( )] ( ) ( ) ( ) ( [ 2 1 ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( 2 1 ) (2 1 2 1 2 1 2 2 2 2 1 1 1 1 2 1s s r r r r s r s r s r s r x xa b b a b a b a (2 22) The electron density for the system can be constructed using Equation 221 with integration over both spatial and spin components: 2 1 2 1 2 2 2 1 2 2 2 1 2 1 2 1 2 1) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( 2 ) ( r r s s dr r r r r dx ds rb a a b b a (2 23) Equation 223 shows the prefactor of 21/2 is normalization constant. In a similar manner, the wavefunction for a system of two electrons with opposite spins occupying the same orbital is given as: )] ( ) ( ) ( ) ( )[ ( ) ( 2 1 ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( 2 1 ) (1 2 2 1 2 1 2 2 2 2 1 1 1 1 2 1s s s s r r s r s r s r s r x xa a a a a a (2 24) For more than two electrons, the determinant is expanded w ith an additional column for each new single particle wavefunction and an additional row for each new electron:
38 ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( 1 ) ( 1 2 / 1 2 / 0 0 1 1 1 2 / 1 1 1 2 / 1 1 0 1 1 0 1 N n N N n N N N N N N N ns r s r s r s r s r s r s r s r N x x (2 25) The notation gets a bit cumbersome at this point. The nomenclature and short hand used by Finnis helps matters.21 One convenient condition for illus trative purposes is to only consider closed shell system where each orbital is doubly occupied. The Energy for the system is given as: 1 2 / 02N n nE (2 26) The prefactor of 2 is the occupanc y. Similarly, the electron density for the many elect ron wavefunction is given by: 1 2 / 0 22 ) (N n nr (2 27) The Slater determinant introduces an electron correlation energy called exchange correlation, where the Coulomb interaction between electrons with parallel spins is reduced compared to the repulsion between electrons with opposite spins.64 This can be illustrated by solving the Coulomb integral between two electrons described by the wavefunction in Equations 222 and 224. In the case where electrons have the same spins, the repulsion between them is:
39 ab ab b a b a b a b a b a b a b aK J dr dr r r r r r r dr dr r r r r dx dx s r r r s r dx dx s r r r r s r r dx dx s r r r s r dx dx r r 2 1 2 2 2 1 1 1 2 1 2 2 2 1 2 1 2 1 2 1 2 1 2 1 2 2 2 2 2 1 2 2 2 2 2 1 2 1 1 1 2 1 2 2 2 2 2 1 2 1 2 1 2 1 2 1) ( ) ( 1 ) ( ) ( ) ( 1 ) ( ] ) ( ) ( 1 ) ( ) ( ) ( ) ( ) ( 1 ) ( ) ( ) ( 2 ) ( ) ( 1 ) ( ) ( [ 2 1 1 2 1 (2 28) Jab represents the Coulomb integral between single particle states and, Kab the exchange integral between states. The same calculation for two electrons with opposite spins occupying the same spatial state: a a a a a a aJ dx dx s r r r s r dx dx s s r r r s s r dx dx s r r r s r dx dx r r ] ) ( ) ( 1 ) ( ) ( ) ( ) ( ) ( 1 ) ( ) ( ) ( 2 ) ( ) ( 1 ) ( ) ( [ 2 1 1 2 12 1 2 1 2 1 2 1 2 2 2 2 2 1 2 2 2 2 2 1 1 1 2 1 2 1 2 2 2 2 2 1 2 1 2 1 2 1 2 1 (2 29) Comparing the two results shows how the Coulombic interaction between electrons is reduced when spins are parallel. The exchange correlation only affects the interaction between electrons with the same spin. 2.2.5 The HartreeFock (HF) method The first useful method for approximating the solution to the many electron Schr dinger equation based on the variational principle was developed by Hartree.65
40 The Hartree method does not solve the true Schr dinger equation but instead solves the equation for noninteracting electrons. In the Hartree approximati on, the electronelectron interaction is approximated with an external potential that represent the mean field of the electron density of the system. The external potential augments the electronnuclear potential to give an effective potential, Veff, acti ng on the electrons. This mean field approximation simplifies the many body calculation in two ways. First the many body Hamiltonian can be treated as a sum of single particle Hamiltonians, Hi: E HN i i 1. (2 30) Secondly, the single partic le Hamiltonian is reduced to a potential and kinetic energy term: ) ( i eff i ir V T H (2 31) Veff contains both the external potential from the nuclear point charges and the potential due to the mean field of the other electrons: ) ( ) ( ) ( dr r r r r V r Vext eff (2 32) where, (r) is the electron density due to all other electrons. The second integral on the right of Equation 232 is defined as the Hartree potential VH(r). Since VH(r) depends on the existing electron density; the Hartree Hamiltonian must operate on a wave function upon which it depends. This means that the problem must be worked iteratively, with successive iterations giving improved wavefunctions, which in turn, can be used in the next iteration. If all goes well, both the wavefunction and the energy conver ge to invariant points where the difference in energy between iterations is suitably small that the solution can be considered converged.
41 The Hartree method is outdated but introduces a few significant innovations that are still indicative of quantum cal culations. The method introduces an iterative solution based on the variational principle. This same technique remains in use for all methods that are solved variationally, which constitutes the vast majority of first principles calculations. The iterati ve nature of the calculations is also the root of the heavy computational expense of first principles methods. The Hartree method also makes use of the mean field approximation which has also carried on in various forms in todays methods. The main critic ism of the method is the Hartree energy does not include the effects of electron correlation. We can at least address exchange correlation by using Slater determinant wavefunctions. This model is called the Hartree Fock (HF)66 method for which the Hamiltonian, F is modified to include the exchange interactions: G V T Fext (2 33) The operator G contains all the electronelectron interactions and is defined as: K J G 2 (2 34) J and K are the Coulomb and exchange operators. The final HartreeFock energy, EHF is found by operating F on a Slater determinant wavefunction: x H eff N n HFE dr r V r dr r V r n n E ) ( ) ( 2 1 ) ( ) ( 2 1 22 1 2 / 0 (2 35) The second term on the right of Equation 235 is the energy of interaction between each electron and the effective core potential defined in Equation 232. The third term corrects for the double counting of the interaction between each electron wit h the Hartree potential that is included in Equation 232. The first term in Equation 235 is the kinetic energy of each single particle state written in Dirac notation.
42 The Dirac style of notation simplifies the expression of many electron integrals.53 The style is also called bracket notation where the ket, n represent the single particle state n. The scalar product of n with the bra n results in the complex conjugate of the single particle state. One other useful scalar product is between n and a vector in real space r : ) ( r n rn (2 36) The last term in Equation 235 is the exchange energy defined as: '1 2 / 0 1 2 /drdr r r n r r m m r r n EN n N n m x (2 37) The integration in Equation 237 runs over four indices, which means the method scales as N4 for N basis functions. This limiting factor restricts the method to small systems. Still, EHF captures the essential qualities of chemical bonding. The main deficiency is that only exchange correlation is explicitly included in the result. However, post HF methods can correct for other correlation effects with chemical accuracy in small systems. 2.2.6 Key aspect s of chemical bonding in molecular orbital theory Wave function theory (WFT) gives the energy of a system of electrons based on molecular wavefunctions, also c alled molecular orbitals, which, in turn, are formed from a linear combination of local single particle orbitals. This molecular orbital based view is referred to as molecular orbital theory. Many of the properties we associate with chemical bonding: orbit al hybridization, bond orders, partial charges and polarization to name a few, can be explained from the perspective of WFT in the context of molecular
43 orbital theory. A few of these key aspects of chemical bonding have to be refined a bit further to subs tantiate the approximations made in the empirical COMB formalism. Since this work deals mainly with the development of variable charge potentials, the first useful concept to explore is the idea of partial atomic charge. WFT provides a means of determini ng the partial charge of each atom based on the degree to which the local atomic basis functions contribute to the molecular wavefunction. In Dirac notation, a single particle orbital can be described as I where the index I refers to the atom and refers to a single particle orbital on that atom. As shown previously, a single particle state can be expanded as a linear combination of basis functions, in which case, expansion coefficients follow the same indexing: I n II C n (2 38) The total number of electrons in the system can be summed up in terms of the expansion coefficients as: electrons n I J I J I J I n J n I n n I nS C C f C f N0 2 (2 39) In Equation 229, fn is the number of electrons in each orbital ranging for zero to two. S is the overlap integral between basis functions: J I SJ I (2 40) The notation in equations such as Equation 229 can get complicated. One useful concept is the density matrix30 whose elements are defined as: occupied n n J n I n J IC C f0 (2 41) Using the density matrices, Equation 239 can be r ewritten as:
44 electrons n I J I J I J I J I I IS N0 (2 42) The charge of each atom depends arbitrarily on how the contributions of the second summation in Equation 242 are divided between atoms. The simplest approach was proposed by Mulliken, who divided the second term evenly between all of the atoms upon whi ch the basis functions resides.6770 The Mulliken charge at each site is then: J J I J I I I IS q (2 43) The methods runs into problems when the orbitals not orthonormal. Often the basis functions on individual atoms are constrained to be orthonormal. However, enforcing the same constraint on the set of molecular orbitals for the system is not mathematically convenient and often not beneficial. In such a case, the summations in Equation 243 do not equal N. Several schemes to address this concern have been proposed, all of which transform the orbitals into an orthonormal basis set. How this transformation is performed differentiates methods such as the L wdin population analysis71 and the natural population analysis (NPA) .72 While these methods maintain charge conservation, the various processes of rendering an orthonormal basis set alters the a bsolute value of charge. Consequently, results of various methods are not comparable to one another. The Mulliken population analysis and its modifications have several other critical deficiencies, which restricts their overall general utility. First, as is apparent in Equation 2 43, the charges depend upon the basis functions so that charges calculated with different basis sets are not comparable to on another. Consequently, the analysis does not provide an absolute value for charge, but instead yields a relative charge subject to
45 the conditions of the calculation. Secondly, splitting contributions evenly between atoms is a rather crude approximation. One could imagine in a real system, with differences in electronegativity between atoms, that the distr ibution is more complex. For these reasons, partial charges are not useful by themselves as quantitative charge values within a material. Their best use is to compare changes in the relative partial charges in response to some change in the system. Occasi onally, a variable charge potential will be touted for reproducing the partial charge generated by a Mulliken or similar population analysis. Such claims are somewhat specious, since the actual value of the charge can vary with the basis set used, and ,while, relative trends in the population of orthogonalized, single particle basis functions may provide qualitative insight, there is no way to ensure relative values quantitatively represent change in the actual charge distribution in real systems. Another concept that is useful for this work and can be derived from the density matrix is the bond order ( BO ), which is a measure of the number of chemical bond that exists between a pair of atoms.73 The range of values begins at zero when no bond exists, and increases by one for each additional bond formed between atoms. Bond strength generally increases with bond order. For example the carboncarbon triple bond has a bond enthalpy of 8.51eV and a bond order of three versus 3.69 eV for the C C single bond with a bond order of one. Hence, it is a critical concept to model when crafting a reactive analytical potential, especially for hydrocarbons where materials properties and chemistry are greatly influenced by the bond order. From WFT, the bond order is simply a sum of the off diagonal matrix elements of the density matrix : occupied J I n n J n I nC C f BO, 0 (2 44)
46 An example used by Finnis to illustrate this concept is the bonding between diatomic hydrogen.21 Each H atom contains one electron in one s orbital, A and B, which combine to form a bonding molec ular orbital, a and an antibonding orbital, b: ) ( 2 1 ) ( 2 1B A b B A a (2 45) In the ground state, a is occupied with 2 electrons, so f0=2 and b is vacant, f1=0. The coefficients for the bonding state,0 AC and0 BC are equal to 21/2. For the antibonding state, the coefficients are 2 / 1 12AC and2 / 1 12 BC. Putting these values into Equation 2 44 gives a BO =1 for the ground state meaning a single bond exits between atoms in the ground stat e. When one electron is promoted into b so that f0=1 and f1=1, then BO =0 and no bond exists. This example also illustrates the essence of chemical bonding, where as atoms interact, the basis functions on each atom combine to form molecular orbitals. The molecular orbitals are both higher and lower in energy than the atomic states. The orbitals are populated from lowest energy states first. If more electrons are found in states lower in energy than the atomic states, then it is energetically favorable to form a bond. The splitting of states continues as more atoms interact until, in condensed matter, the states broaden into bands.74 Within a crystal, if we consider the exclusion principle, each band constitutes a single particle spin. The bands represent the distribution of energy states within the solid phase which can be described using the density of states
47 (DOS). If the single particle states, n and energies, n, are known, the global DOS for the system, D() is defined as: n nD ) ( ) ( (2 46) whe r e, is the delta function. The function D() represents the number of states with energies in the range of +d The band energy, EBand, is the sum of t he energies of the states given by the integral: n E n n BandFd D f E ) ( 2 (2 47) At zero temperature, the upper limit of the integral is the Fermi energy, EF. The factor of two accounts for a closed shell system with each state being doubly occupied elect ron of opposite spin. A local DOS for each orbital on any atom can be defined as a function that measures the relative contribution of each state of energy, to the total DOS of the orbital. This is achieved by first defining the relative weight of each of molecular orbital: ) ( ) (n n n J n I J IC C D (2 48) The global DOS is then the sum if: J I I J J IS D D ) ( ) (. (2 49) The local DOS for orbital I can be pulled from Equation 249 as: J I J J I IS D D ) ( ) (. (2 50)
48 With an orthogonal basis, DI() is simply DII(), the sum of which over gives the DI() for each atom. The local DOS provides a route to approximating bond energies in condensed matter, which will be used to derive the short range contributions in COMB. 2.3 Density Functional Theory (DFT) of Inter Atomic Interactions The main limitation with the single particle HF methods is the error in the correlation energy. This error can be corrected with any of several post HF methods. The post HF corrections are naturally expensive and scale rather poorly with the number of electrons, N4 for the mini mally useful methods .58 Unfortunately, the accuracy of the methods only improves in direct proportion with the expense to the point where the most trusted post HF methods, coupledcluster theory, scales as N8 where N is the number of basis functions. The poor scalability restr icts the calculations for fairly small systems containing at the most a few dozen atoms. Typical calculations are performed on systems on the order of one to ten atoms. One solution to the many body correlation problem was advanced by Thomas and Fermi who found that for a system containing N electrons, there exist a unique functional of the electron density E[r], which when minimized with respect r gives the true ground state energy of the system.75, 76 A functional is a function that takes a function as its arguments. In this theory, the energy f unctional takes the electron density as its argument hence the name density functional theory. The appeal of this theory is the many body correlation effects are inherently included in the electron density. The resulting ground state energy is the exact ground state energy that can be solved as long as the appropriate functional is known. Unfortunately, the exact functional is not known and the theory in its current state relies on empirically derived
49 approximate functionals, thus yielding an approximate solution to the exact ground state energy. Still DFT includes correlation effects in a comparatively efficient manner, so can be applied to larger systems with reasonable cost. Systems of several hundred atoms are tractable on current computer systems. This has lead to the acceptance DFT as the predominant first principles method in computational materials science, especially when electron rich atoms such as metals are involved in the calculation. 2. 3.1 The KohnSham Functional Modern empirical potenti als for materials stem primarily from DFT and follow an evolutionary development along the lines of tight binding theory (TB DFT). Tight binding theory is semi empirical density functional theory, where approximations are made to the total energy functional to increase the efficiency of the calculation. The theory encompasses of a wide range successively sever approximations, which result in correspondingly faster calculations, but at the cost of progressively limited transferability. Extrapolating this tr end to its limit leads to fully classical empirical potentials. The most useful tight binding energy functionals for deriving empirical potentials are themselves based on the KohnSham functional which can be broken down into its component energy functionals as:27, 28 ZZ eZ xc H s KSE E E E T E (2 51) The main challenge with energy functionals preceding the KohnSham approach was in dealing with the f unctional derivative of the kinetic energy. Previous methods approximated the kinetic energy, which proved to be a significant source of error. The KohnSham approach works around the approximations to the kinetic by introducing a
50 reference system of noni nteracting electrons with a density of The kinetic energy for this reference system is given by the functional, Ts which is the first term in Equation 2 51. Ts does not yield the true kinetic energy, but rather the kinetic energy of a reference system. The actual kinetic energy funct ional requires a correction term, T=Ts+ T, which is assumed to be included in the Exc. The other functionals give, from left to right, the Hartree energy, the exchange and correlation energy, the electronnuclear interaction and the nuclear nuc lear interaction. The annotation is modified for simplicity where (r) is abbreviated as just which will be used when appropriate from here on. The ground state for any energy functional is determined by finding the minimum in the energy functional w ith respect to electron density. In a manner similar to functions, the critical points for a functional correspond to point where the functional derivative is equal to zero. The ground state density is, therefore that which corresponds to the point where the functional derivative of EKS with respect to density is zero with the constraint that total charge in the system is conserved: ) ( ) ( ) ( r V r V r V T Eext xc H s KS (2 52) is the Lagrange multiplier used to enforce the constraint of charge conservation meaning the following condition is satisfied: N dr r ) ( (2 53) is an important concept that essentially represents the change in energy with electron densi ty or in terms of charge, the change in energy as an electron is added or removed
51 from the system. The other potentials in Equation 252: VH(r), Vxc(r) and Vext(r), are the functional derivates of the corresponding energy functionals. Solving the KohnSh am equation requires a bit of clever manipulation. The solution requires defining a reference system of non interacting electrons as: ) ( r V Teff s (2 54) When Veff, the effective potential acting on the electrons is defined as: ) ( ) ( ) ( ) ( r V r V r V r Vext xc H eff (2 55) Equations 252 and 254 are equivalent and are solved with the same ground state density and kinetic energy. However, by choosing a system of noninteracting electrons as the reference system the, exchange and correlation contributions can be dropped and Veff becomes essentially an external potential, Vext. This sets up an eigenvalue problem which can be solved variationally: ) ( ) ( ) ( 2 12r r r VKS n KS eff (2 56) Eigenvalue equations in this context are referred to as the KohnSham equations The wavefunctions, ks are the single particle KohnSham orbitals that correspond to the ground state solution of the Kohn Sham equations. The electron density is given by: 2) ( ) (n n nr f r (2 57) The summation runs over all electrons. The term fn is an integer that represent the orbital population ranging from zero to two. The main advantages of this solution are that the ground state density and kinetic energy are determined via a self consistent solution to a single particle eigenvalue problem. The problem a sequence of iteration
52 since Veff is dependent on Most other approximations are gathered together in the Exc. The method does not side step the challenges with the exchange and correlation functional but rather lumps all approximations together into one empirically parameterized functional. Finnis reworked EKS into a form that is useful for deriving empirical potentials.21 First, the potentials such as Veff can be interpreted as operators effV giving the following expressions: n n sn T n f T (2 58) n eff n effn V n f dr r V r ) ( ) ( (2 59) Combining the two expression gives: dr r V r n V T n f Teff n eff n s) ( ) ( (2 60) Replacing Veff with its comp onents gives Finniss form of the functional: ZZ XC XC H n eff n KSE E dr r V dr r V n V T n f E ] [ ) ( ) ( 2 1 ] [ (2 61) The second term in the equation corrects for double counting of the Hartree energy since it is counted twice in the Hamiltonian. The third term subtracts any exchange and correlation contributions from the Hamiltonian and includes them in Exc. 2. 3.2 The HellmanFeynman theorem and inter atomic forces If we intend to use the potential for dynamic simulations or optimizations then we have to deal with forces acting on the atoms, w hich in quantum mechanics, are found through the HellmannFeynman theorem.77 The principle of the theorem states that
53 change in the energy due to change in some factor that alters the exter nal potential is given by: 0 0 R V R Eext (2 62) The change in energy is not dependent upon the change in electron density at all; only the change in the external potential. As an example, provided by Finnis ,21 consider a system in the ground state where the energy is given by s ome functional, E. If the external potential is alt ered with some small perturbation, such as a small change in interatomic displacements the resulting change in energy to first order is: ) ( ) ( ) ( ) ( ) ( )] ( [ )] ( [ r dE dr r dV r dr r V r d dr d d r dF r dEZZ ext ext (2 63) For clarity, F represent s a ll other components of the functional derivative. With the adiabatic approximation, the electron density is minimized to the ground state for each small perturbation in the Vext. Since ris minimized, all terms that depend on dr add up to zero leav ing: ) ( ) ( ) ( )] ( [ r dE dr r dV r r dEZZ ext (2 64) Equation 264 restates the HellmannFeynman theorem, where the change in energy is not dependent upon change in the electron density. ) ( ) ( ) ( )] ( [ r dV dE r r dV r dEext ZZ ext (2 65) This theorem greatly simplifies the task of calculat ing inter atomic forces. Once the ground state density is known, the forces acting on the nuclei in response to the ground state charge distribution can be solved through classical electrostatics. The self
54 consistent quantum calculation is only necessary to find the ground state electron density distribution. 2. 3.3 The second order Kohn Sham functional T he HellmannFeynman theorem can also be turned around to derive functionals that yield the change in energy and electron density in response to a change in Vext by introducing a scalar parameter, that determines the external potential over the range of the change between Vext1 and Vext2. 21, 78 ranges from zero to one. The potential in terms of is: ) (1 2 est est est estV V V V (2 66) The HellmannFeynman equation that corresponds to this condition is : ZZ ext extdE dr d r V r V r r dE ) ( ) ( ) ( )] ( [1 2. (2 67) Integration with respect to gives the energy functional. However, solving the integral requires an approximation for the dependence of (r) on The first order approximation is to disregard any dependence (r) on 78 which yields the following first order functional: 1 2 1 2 1 1 1 2 ) 1 () ( ) ( ) ( ] ), ( [ ] ), ( [ZZ ZZ ext ext ext extE E dr r V r V r V r E V r E (2 68) The error in the approximation varies linearly with Vext and so is only reliable for small variations in Vext. However, the functional does not depend on a change in electron density, and, so, can be solv ed efficiently in situations where it is applicable. An approximation to second order can be made by assuming the dependence of (r) on is linear.21, 78 The resulting functional is essentially obtained by replacing 1(r) in the first order functional with the mean density [1(r)+2(r)]:
55 dr r V r V r r E E dr r V r V r V r E V r Eext ext ZZ ZZ ext ext ext ext) ( ) ( ) ( ) ( 2 1 ) ( ) ( ) ( ] ), ( [ ] ), ( [1 2 2 2 1 2 1 2 1 1 1 2 ) 2 ( (2 69) The KohnSham analog to the second order functional is given as: )' ( 2 1 )] ( ) ( [ )] ( [ )] ( [2 1 ) 2 ( r r C dr r V r V r E n H n r Ein in Hin XC in in in XC in in (2 70) The superscript in indicates a reference input densi ty, in. Functions and functionals that are dependent upon in are also tagged with the in superscript. The first term on the right of Equation 270 is the sum of the single particle KohnSham eigenvalues. The premise of the second order functional is that the energy is determined as function of (r) where (r) is in(r) plus a small, continuous change in (r), (r)=in(r)+(r). Likewise, the single particle states are assumed to change smoothly with small changes in (r). The first term is a sum of the eigenvalues found by operating the following Hamiltonian on the single particles states, n: in ex in H ext inV V V T H (2 71) Harris analyzed Equation 270 further focusing on the last term, which is comprised of the second order te rms resulting from the expansion of the EH[(r)] and Exc[(r)] .79 inr E r r r r Cxc in )] ( [ 1 ) (2 (2 72) Solving the second order functional requires solving a KohnSham equation of the form n n V Tn eff ) () 2 ( (2 73) The effective potential i n this case is :
56 ) ( ) ( 2 1 ) ( ) ( ) () 2 (dr r r r C r V r V r V Vin in XC in H ext eff. (2 74) Equation 274 indicates that the quantity defined in Equation 272 represents the change in the Hartree and exchange and correlation potentials due to change in the electron density. Equation 270 is a convenient functional from which to begin the derivation of a classical potential that gives the energy as a function of change in the nuclear positions. The functional requires a self consistent step to arrive at the ground state density. However the error in the energy is shown to be third order so small variations around the true variational solution are tolerable. 2.3.4 The first order KohnSham functional Harris, Foulkes and Haydock found a first order approximation of the second order KohnSham functional can be made by basically dropping the second derivative term from Equation 270:78, 79 ] [ )] ( ) ( [ ] [ ] [2 1 in ZZ in H in XC in in in XC n in HFE dr r V r V E E (2 75) The appeal of this functional is that it is not variational. The eigenvalues are determined once based on the input density and so determining the ener gy only requires one step to solve the KohnSham equation. Furthermore, the error in energy is second order in terms of the difference between the input and exact charge densities. With a careful choice of in and Vxc, the functional can actually be more accurate than the full self consistent solution. The second order error is satisfactory for many applications considering the efficiency gained by avoiding a self consistent iterative procedure. The consequence of not being variational is the calculated energy may be higher or lower
57 than the actual ground state energy. There is no lower bound. Hence the results are only as reliable their validation allows. 2. 4 A Self Consisten t Classical Electrostatic Pote ntial The purpose of an analytical potential is to determine the energy and forces acting on the nuclei as a function of the position of the nuclei. The first and second order functionals can serve as a starting point to derive an empirical potential with a foundation in density functional theory. The first order Harris Foulkes functional provides a theoretical basis for several fixed charge, many body potentials. Brenner has shown how the Abel Tersoff and REBO potentials have their theoretical origins in this functional.52 Likewise, Finnis has derived Pettifors analytical bond order potential ,21, 42 the Finnis Sinclair potential21 and a basic ionic potential starting with the Harris Foulkes functional.21 The second order functional provides a basis for a self consistent tight binding model, in which several electrostati c analytical potentials have their theoretical basis. In particular, models that include change in the electrostatic interactions in response to the local environment such as a polarizable core shell models ,80 can be derived from the second order functional. The dynamic charge potential central to this work also has its origins in the second order functional as well. If the input charge density in second order Kohn Sham functional is equal to the exact Kohn Sham ground state charge density, the functional may be rewritten as: ) ( 2 1 )] ( [ )] ( [) 2 ( r r C r E r Ein HF in. (2 7 6 ) The first order portion of the functional is equivalent to EHF, the first order Harris Foulkes functional. Equation 27 6 breaks the task of calculating the energy into two par ts. The
58 Harris Foulkes functional gives a nonself consistent solution to the energy that is only dependent upon the input charge density. The last term on the right of Equation 230 gives the change in the Hartree and exchange and correlation energies to second order. Together both terms provide a basis for a reactive analytical potential. Equation 27 6 also allows the derivation to be separated into self consistent and nonself consistent components of the potential. 2.4.1 Self consistent electrostatic energy contributions Ideally, a density function that can capture the change in energy with respect to electron density with one atom specific variable would solve the problem in as efficient a manner as possible. The partial charge on each atom can repr esent the change in density as long as the density distribution is rigid, meaning its volume and shape are fixed. Such a function is provide by Streitz and Mintmire in their electrostatic plus model (ES+):22 ) ( ) ( ) ( ) ; (i i i i i i i ir r f q r r q r (2 7 7 ) is an effective core charge ranging between 0 and the total number of valence electrons in the atom. f(r) is a function that describes the radial charge distribution. If we also neglect the change in Exc and any covalent ef fects, then the change in energy as a function of charge is a sum of the change in energy in each atom. With this rigid ion approximation, the change in energy becomes the sum of the Coulombic interactions between charge densities. ) ( ) ( 2 1 )] ( [ drdr r r r r r Ees (2 78)
59 With rigid charge density distributions, Equation 27 8 can be thought of as the classical electrostatic interaction between charge densities the sum of which gives the total electrostatic energy, Ees, of the system. Solving the Coulombic integr al with the charge density function of Equation 278 yields several electrostatic potentials. The columbic interaction between charge densities, Eqq is given as the integral over charge density distribution functions, f(r): j i j i j j i i qq ij i j j qq ij i qqdr dr r r r r f r f J q J q E ) ( ) ( 8 10 (2 79) 0 is the vacuum permittivity. The energy contribution due to the electronnuclear interactions, Eq, is determined as the Coulombic attraction integrals between the valence charge and the effective core charge: j i j i j j i i i i i i j i ij q ij i j j qZ ij i qZdr dr r r r r f r f dr r r r f q r J Z J q E ) ( ) ( ) ( ) ; ( 4 10 (2 80) The energy due to nuclear nuclear interactions is given by: r dr r r r f dr r r r f dr dr r r r r f r f r J Z J Z Ej j j j i i i i j i j i j j i i j i ij ij i i j j ZZ ij i ZZ1 ) ( ) ( ) ( ) ( ) ; ( 8 10 (2 81) When the core charge is approximated as a point charge and the densities are rigid, EZZ is a constant with respect to charge. Therefore, we can assume EZZ is included in the cor e core interaction of the Harris Foulkes energy.
60 Lastly, since the summations run over all atoms, there are intraatomic Coulombic interactions between the valence and the core on the same atom,q iiJand the electronelectron interaction on the same atom, qq iiJ These onecentered integrals comprise the electrostatic self energy. With the rigid atom approximation, the self energy becomes a sum of atomic self energies, Vself: i i qq ii i i qZ ii i self i selfq J q q J V E 2 1 4 10 (2 82) Incorporating Equations 280 to 282 in Equation 276 gives a second order functional as: qZ qq self in HF inE E E r E Q r E )] ( [ ] ), ( [) 2 ( (2 84) The choice of charge distribution function is somewhat arbitrary. Ideally, the function should be computationally efficient yet reliably reproduce properties of the material. The simplest approximation would be to use a point charge for the density distribution, which is common in potentials for ionic materials when efficiency is paramount. However, the point charge approximation has significant erro r when atoms are separated by close distances and goes to infinite values as atoms approach too closely together. This Coulomb catastrophe is avoided with spherical charge densities, which, instead, go to finite values as densities overlap. The effect of using spherical density functions as opposed to point charges is illustrated in Figure 2 1 which plots the interaction between two charged atoms versus separation distance. In Figure 2 1, dashed lines represent a point charge model and solid lines represent Coulomb integrals over spherical charge distributions. The charge on each atom is fixed at 1 with like charges giving positive values and opposite charges giving negative values.
61 The functions are nearly equivalent at long separation distances, but deviate near the typical range for bond lengths in solids. The point charge model goes to infinity at close separation which leads to a catastrophic failure of the model. Spherical distribution functions have been shown to be effective in tight binding models and are the mathematically convenient.52 A rigid spherical charge density simplifies the electrostatic portion of the Equation 276 Spherical densities also simplify evaluation of the Harris Foulkes energy, which will be addressed in the next section. Consider ing these advantages, spherical densities strike the best balance between efficiency, reliability and performance. Streitz and Mintmires spherical density function is mathematically convenient and has been shown to perform well in analytical potentials.22, 40 The function has one atom specific variable, which controls the radial rate of decay of the charge density: ) 2 exp( ) (1 3 i i i ir r r r f (2 84) Incorporating Streitz and Mintmires density distribution function into the interaction integrals yields the following exponentially decaying functions. When i=j: 2 3 26 1 4 3 8 11 1 ) 2 exp( 1 ) (ij ij ij ij ij qq ijr r r r r r J (2 85) a nd, when ij, 3 3 4 2 3 3 4 2 2 2 4 2 2 4) 2 exp( 3 ) 2 exp( 3 ) 2 exp( ) 2 exp( 1 ) (i j i j ij ij j i j j i j i ij ij i j i i j i j ij j i j j i j i ij i j i ij qq ijr r r r r r r r J (2 86) The charge nuclear interaction integral yields:
62 ) 2 exp( 1 ) 2 exp( 1 ) (ij j ij ij j j ij i i i ir r r r dr r r r f (2 87) The charge charge interaction, Jqq has an r1 term. Summations over r1 functions are conditionally convergent, which leads to problems in calculating the total energy. A few techniques have been developed to deal with such summations. For this work, the sum is solved via a direct Wolf summation over charge neutral spheres.81 The direct Wolf summation is computationally efficient and is easily incorporated into parallel codes. The technique adds an additional q2 dependent self energy term. This additional atomic self energy is assumed to be included i n the self Coulomb interaction, qq ijJ ; otherwise, the self energy is unnecessarily dependent upon parameterization of the long range summation. The self energy terms in the rigid atom model comprise the energy to form a charge on each atom, whic h as shown by Rappe and Goddard,14 may be approximated through a Taylor series expansion of the energy as a function of charge: 2 2 2 02 1 ) (i i i i i i Self iq q E q q E E q V (2 88) Truncating at second order and solving at neutral charge as well as q=1 and q= 1 gives the following relationships: 0) 0 (i Self iE V (2 89) 2 2 02 1 ) 1 (i i i Self iq E q E E V (2 90) 2 2 02 1 ) 1 (i i i Self iq E q E E V (2 91)
63 VSelf(1) is defined as the first ionization potential (IP) and VSelf( 1) is the negative of the electron affinity (EA). From the above equation we can get: i iEA IP q E ) ( 2 1 (2 92) qq ii iJ EA IP q E ) (2 2, (2 93) The first derivative of energy with respect to charge is the equivalent to the Mulliken electronegativity, With the rigid atom approximation, is equivalent to one centered electronnuclear interaction, JqZ. The second derivative of energy with respect to charge is Jii, the Coulombic interaction between valence electrons i n the same atom. Equations 2 92 a n d 293 establish a way of determining Vself directly from calculated or experimentally measured ionization potentials and electron affinities. When the expansion is truncated at second order, the resulting function is quadratic and exactly represents the first ionization potential and electron affinity. For many elements, the higher ionization potentials vary significantly from the quadratic solution. Additional higher order terms may be added to the function to fit higher oxidation states. In Equation 2 82, where Vsel f is defined as the one centered integrals, the value is divided the vacuum permittivity (40)1. When the coefficients are derived from empirical val ues, such as with Equations 292 and 293, the permittivity is assumed to be included in the measured value and is not included explicitly in the equation. Equation 288 provides a classical approximation to the change in energy with respect to the partial charge on each atom. The self consistent problem is reduced to determining the partial charge on each atom that corresponds to the ground state electron density configuration.
64 2.4.2 Dynamic charge equilibration As defined previously, the electronegativity is the partial derivative of energy with respect to charge ( E/ qi)qj. In Equation 25 2, the partial derivative of energy with respect to electron density is set equal to the Lagrangian multiplier that enforces the constraint of charge conservation. With the rigid atom approximation, where variation in electron density is only a function of charge, the derivative with respect to electron density can be truncated to the derivative with respect to partial charge. Likewise, it has been shown by Parr, that the electronegativity is also equivalent E/ (r) when the (r) is the self consistent KohnSham density .30 By Equation 252, the ground state corresponds to the condition where the electronegativities are equal to at all sites, provided that the total charge in the system is conserved. Likewise, the distribution of partial charges that corresponds to the ground state electron density distribution satisfies the condition that electronegativity is equal at all sites. The partial charges are selected as those that satisfy this unique electronegativity equilibration condition (EE) with the constraint that total c harge is conserved. The most computationally effective means of enforcing the constraint is though an extended Lagrangian method similar to that of Rick et al. ,82 where the constraint of charge conservation is enforced with an undetermined multiplier. The Lagrangian for the system is: N i i i qZ qq self in HF N i i Qi N i i iQ E E E r E Q M r m L1 1 2 1 2)] ( [ 2 1 2 1 (2 9 4 ) Here, Mq is a fictitious mass of each charge, m i s the mass of each atom, is the undetermined multiplier.
65 From the Lagrangian, we can derive equations of motion for both the real space coordinates and the dynamic charges as: i i ir r Q E r m ) ( (2 95) i i i QQ r Q E Q M) ( (2 96) Since the total charge is conserved and therefore a constant motion, it follows that: i iQ 0 (2 97) Solving the above two equations shows to be the negative mean electronegativity for the system: i iN 1 (2 98) This leads to an equation of motion in charge space as: i i QQ M (2 99) Equation 299 sets up a scheme to determine the partial charge on each atom concurrently with the spatial coordinates during the course of a dynamic simulation. The charge equations of m otion and be solved following each real space projection of the atomic positions using standard numerical integration algorithms. For this work, a velocity Verlet algorithm83 is used to integrated both the real and charge space equations of motion. In this scheme, the total energy of system is not conserved. The total force on each atom can be determined as the negative total derivative of the potential at each si te. As an example, the component of force in the x direction acting one any atom is given by:
66 i i i xx E x q q E dx dE f (2 100) When electronegativity equilibration condition is exactly satisfied so that E/qi= at all sites and the total charge on the system is equal to zero, the first term on the right of Equation 2100 is equal to zero since: 0 i j i iq x r q q E (2 101 ) Therefore, if the electronegativity and charge conservation conditions are rigidly sa tisfied, we do not have to determine derivatives of charge with respect to position. In the scheme presented here, the EE condition is never exactly satisfied. The exact solution requires, first, a linear derivative of the energy with respect to charge and, secondly, solving N linear equations for the charge on each atom. The formalism presented here does not have a linear first derivative. With dynamic charges, fluctuation about the time average leads to a nonconserved condition at each time step. Although absolute energy conservation is not achievable with such a scheme, conservation within reasonable limits is possible by iteratively projecting the charges with damped velocities at each time step until the electronegativity at each site is nearly equal w ithin a determined tolerance: i i i Qq Q M (2 102) is an empirically optimized damping factor The solution in Equation 2101 only applies to charge neutral systems. Its a convenient to maintain the constraint of charge neutrality for that reason and also because it avoids discontinuities in the energy when determined over periodic cells.
67 2. 4.3 Analytical b ond order potentials for multiphase simulations Remaining part of the basic potential formalism requires an analytical approximation for the Harris Foulkes energy. As noted earlier, the analytical bond order potential formalism has been shown to adeptly apply to many different bonding environments. The transferability of that potential comes in part from it origins in the Harris Foulkes functional as shown by Brenner.52 A review of Brenners derivation highligh ts the approximations taken to achieve a purely analytical potential that emulates the underlying quantum physic of chemical bonding. Abell derived a general form of an empirical potential to describe the attractive electronic energy contribution to inter atomic interactions as a sum of bond energies:23 i j ij A ij el ir V b E ) ( (2 103) The function describes short range bond interactions and so is summed over nearest neighbors. The term bij is an empirical correction term that modifies the bond energy based on bond order. Abell was able to show that the bond order is most dependent upon the local coordination number of the atoms, Z and to first approximation b Z1/2. Higher local coordination yields a lower bond order and hence a weaker bond. To be useful the attractive bond order potential must be c ombined with a short range repulsion term which approximates the short range repulsion between neighboring electron densities as atoms begin to overlap. The choice of functions for short range repulsions is somewhat arbitrary. Historically, functional forms are selected based on their accuracy in reproducing the empirical bond data and mathematical efficiency. Abell employed exponential functions for both the short range attraction and repulsive contribution to give an inter atomic potential:
68 i i cohE E i j ij ij ij ij ij ij ij ir b B r A E ) exp( ) exp( (2 104) A, B, and are fitted parameters. Ecoh is the cohesive energy defined as the energy gained per atom when separate charge neutral atoms are bound into a cohesive solid state: )) ( ( 11crystal E E N ETot N i i coh (2 105) To derive energy per atom as given in Equation 2 105, the Harris Foulkes functional (Equation 2 76) can be simplified by again using rigid, charge neutral, spherically symmetric charge densities centered on each atom as an approximation to the input density. In this case, as shown by Foulkes and Haydock,78 the nuclear nuclear and Hartree energies can be reduced to the sum of the Hartree energies of separate atomic densities, which is a constant, CH, and the a sum of pairwise screening potent ials between each atom: i i j ij ij Hr V C ) (2 1 (2 106) Equation 2107 should be corrected for non pair additive contributions due to the exchange and correlation functional, Vnp. However, Foulkes and Haydock78 demonstrate that in regions of minimal overlap between multiple atoms, a pairwise correction term is suitable, in which case Vnp can be added to Vij in Equation 2107. Foulkes and Haydock also show that with spherical charge densities, VXC can be approximated with the expression:78 ) ( ) ( ) ( r U r V r Vi i xc (2 107)
69 V(r) is an additive atomic potential and U(r) is a nonadditive term that corrects for nonlinearity in the exchange and correlation functional. U(r) is a small contribution to the total energy and may be neglected. If U(r) is neglected, the energy for the system given by the Harris Foulkes functional can be approximated as a sum of pairwise interactions plus the sum of single particle eigenvalues based on the input charge density: i i j n n ij totr E ) (. (2 108) As pointed out in Equation 247 in section 2.2.6, the sum of the energies of the single particle states, n, is the band energy, Eband, which is related to the global DOS D() by Equation 247. With the 1st order functional, the band energy can also be defined in terms of Hin as: J I in J I J I bandH E (2 109) The bond energy, Ebond, can be defined in a similar way as the inter atomic contribution to the band energy: I J J I in J I J I bondH E, (2 110) A comparison of Equation 2110 to the definition of bond order given in Equation 244 reveals importance of bond order in determining the bond energies. As with the Eband, the bond energies may also be defined as a function of the local DOS for each single particle state: J J in J I bondFd D H E ) ( ) ( 2. (2 111) Since the objective is an approximation of the energy per bond and knowing that the global DOS can be recovered as a sum of the local densities of state, the
70 appropriate approach is to approximate the local DOS for each eigenstate. DJ() is a distribution function, the shape of which can be characterized by its moments ,84 defined for the pth moment as: Fd DJ p Jv p ) ( (2 112) Comparing equations 2 112 with 2 47 shows the moments related directly to the Hamiltonian through n p n p n pn H n (2 113) and forJv p: J H Jp J p (2 114) The moments are a measure of the shape of the DOS .84 The first moment is t he center of gravity of the distribution given as: FI I Atomic I AtomicD d D 0 ) ( ) ( ) ( ) (1 (2 115 ) With orthogonal orbitals and neglecting charge transfers, the choice of centering the distribution around the energy of the individual atomic orbital, atomic, sets the first moment to zero. The second moment gives the width of the DOS84 defined as: I I AtomicD ) ( ) (2 2 (2 116) The third moment characterizes the degree of skewness in the distribution and the fourth moment relates to the tendenc y to form a gap in the DOS. The cohesive energy of the soli d relative to the free atoms is mainly related to the spread of the orbital energies, which is represented by the second moment. In fact, it can be shown that with approximations, the orbital ener gies on each atom are proportional to the square root
71 of I 2 84 This relationship means the orbital energies can be found without having to determine the higher moments and without having to solve any Schrodinger equations. All that is needed isI 2and the equation of the line that relates 2 1 2 Ito the orbital energies. The advantage of this method is the moments can be determined without explicitly knowing the DOS through the moments theorem ,84 which states that the pth of the local DOS on each atom, DI(), is determined by the sum of all paths between p neighboring atoms that start and end on I For the second moment, the path is just two hops: one to and one back from each neighbor. In other words, the second moment can be determined as the number of nearest neighbors, Z and the orbital energies are proportional to the square root of Z This approximation is referred to as the second moment approximation and forms a basis for several coordination dependent potential energy functions. With perfect crystalline systems, the coordination on each site is clearly defined. The calculation of Z in disordered systems or when defects or interfaces are present gets more ambiguous as clear demarcations between neighbor shells disappear. A flexible definition of what is considered a neighboring atom is required in such circumstances. Finnis and Sinclair found one solution by replacing the sum over neighbors with the sum of function that decays exponenti ally with distance between atoms.21 2 1 I J r Ie E (2 117)
72 is a n empirical parameter. Multiplying Equation 2117 by an empirical proportionality constant, B yields a potential energy function for the bond energy per atom. A short range pairwise repulsive term is needed to balance the bond energy which, when combined with the bond energy gives the Finnis Sinclair potential energy function: I I J I I J r r S F cohe B Ae E2 1 (2 118) A, B, and are empirical parameters. As shown by Brenner, with a little manipulation, Equation 2118 can be recast in a form similar to t h e Abell bond order potential energy function.85 I J J I K r r r Ie Be E2 1 ) ( 21 (2 119) In this form, the bond order, bIJ is: 2 1 ) (1 J I K r r IJe b. (2 120) A significant im p rovement t o the Abell bond order potential was made by Tersoff by including the effects of bond angles and symmetry in the bond order expression.11 i ij i k ik ij IJr r g i b 2 1 ,) ( ) ( 1 (2 121) and are empirical parameters. g() is a function that depends on the bond angle. 2 2 2 2 2) cos ( 1 ) ( h d c d c g. (2 122) is a function that penalizes differences in bond lengths rij and rik. ] ) ( exp[ ) (m ik ij m ik ijr r r r (2 123)
73 c, d, h, and m are empirical parameters. The Tersoff formalism is one of the most effec tive models for modeling covalently bound systems. In COMB, the Tersoff formalism is used as the short range potential energy function. When atoms have zero charge, the bond energy is equivalent to the Tersoff bond energy. Additional modifications are included in COMB for charged states, which will be discussed in the next section. 2.4.4 Beyond the rigid atom approximation The model presented to this point is basically self consistent variable charge electrostatic potential coupled to charge independent short range potential energy function. This approach is equivalent to models such as Streitz and Mintmires ES+ which exhibit deficiencies that hinder their wide spread application. The SM potential, for example, can not model covalent systems and does not predict the ground state energy for alumina, the material it was originally designed to model.86 The limitation of models of this type is due to the rigid ion approximation, which divides the potential energy function into separate electrostatic and bond energy contributions. The bond energy contribution is independent of charge and so is not altered with change in ch arge. Consequently, the bond energy does not change dramatically enough to clearly distinguish ionic, metallic and covalent bonding. At best the bond energy is fit as a compromise between two bonding environments. Such is the case with the SM potential, which is fit to balance metallic and ionic bonding for aluminum. In the COMB formalism, where the goal is to simulate multiphase interfaces, the rigid ion approximation is too severe. The performance of the potential is greatly improved by
74 incorporating a few additional changes in the potential energy due to change in the electron density distribution. 184.108.40.206 Charge dependent short range energy contributions A change in the partial charge on an atom indicates a change in the population of the frontier orbitals. This directly affects the bond order, which is a reflection of the relative population of the frontier bonding and anti bonding states. The change in charge also changes the effective ionic radii, which influence both short range repulsion and the bond energy. A recent potential by Yasukawa addresses theses changes in bond energy with charge.17 The Yasukawa formalism is equivalent to Tersoffs potential at neutral charge, where, the short range energy is composed of a repulsive pair interaction, VR, an attractive pair interaction, VA, and dispersion interactions, Vvdw: The potential deviates from a pure Tersoff potential in that the magnitude and radial decay of the pair interactions vary with charge. Charge dependent correction functions, Di(qi) are added to the exponential decay coefficient of the short range repulsive energy to reflect the change in atomic radius with charge. A similar correction is added to the short range attraction, to reflect the change in overlap with charge: ) ( ) ( exp ) ( ) (2 1 j j j i i i ij ij ij ij c j i ij Rq D q D r A r F q q r V (2 124) ) ( ) ( exp ) ( ) ( ) (2 1 j j j i i i ij ij j i ij ij ij ij c ij Aq D q D r q q B B b r F r V (2 125) In Equations 2124 and 2125, Aij, Bij, ij and ij are Tersoff parameters that are determined for each bond type. The change in short range contributions, Di(qi), is specific to each element type according to: i D i i in i U D U i iq Q b D q D ) ( (2 126 )
75 i i i D n i i iL U U L DQ Q D D b 1 (2 127) ) ln( ln ) ln( lni i i i i i iL U U L U U DQ Q Q D D D n (2 128) DU and DL are parameters that reflect the change in atomic radius with charge. Likewise, QU and QL are the atomic charges that correspond to the limits of the valence shell. The short range attraction contribution is further modified with the function Bij *(qi,qj) that decreases the bond order with increasing charge: 2 1 *) (j i j i ijB B q q B (2 129) i B i i in O i B B iQ q b a B* (2 130) i B BQ a bi B n i i 1, (2 131) 11 i B i in i O BQ Q a (2 132) i iL u iQ Q Q 2 1 (2 133) i i iL u OQ Q Q 2 1 (2 134) The bond order term, bij is equivalent to Tersoff. Also, as in Tersoff, the short range pairwise interactions are terminated by multiplying with a cutoff function, Fc(r) :
76 0 cos 1 2 1 1 ) (S S S S S CR S R r S R R F S S S SS r S r R R r (2 135) In Equation 2135, the interactions are terminated smoothly between two cutoffs, Rs and Ss. With these modifications, the bond energy goes to zero when the charge on each atom corresponds to a filled valence shell. In such cases the attraction between atoms is described entirely by the electrostatic functions, reflecting a purely ionic state. With charges between the valence shell limits, bij contributes a varying degree of covalency to the bond energy, with a maximum contribution at zero charge. 220.127.116.11 Polari zability and the point dipole model Electronic polarization is a distortion in the electron density in response to an external electric field. The inclusion of explicit polarization has proven to be useful in classical simulations where they have been show n, for example, to stabilize the complex crystal structures of some oxides such as alumina. It is suspected that in alumina, the induced dipole and quadruple moments of the anion stabilize corundum lattice.41 Induced dipoles influence the interactions with small polarizable molecules such as water and O2. Con sequently, many current solvent models include explicit polarization effects. The inclusion of explicit electronic polarization is especially pertinent to this work, since two main interests are the study of oxidation and aluminum oxides. The polarization response can be decomposed into several modes. In the classical response, the electron density is distorted from the equilibrium distribution in response
77 to additional potential generated by an external field. This response is generally long ranged and can be approximated with decent fidelity as a sum of atomic polarizations. At short bond distances, two additional modes of polarization comprise a significant portion of the total induced moment. At short separation distances, where the repulsion between electron densities of neighboring atoms is strong, the electron density distribution can be asymmetrically distorted way from the nuclei. This is what Dick and Overhauser labeled the short range interaction mode in their theory of dielectrics.80 The effect increases inversely with separation distance. The second short range mode, what Dick and Overhauser referred to as exchange charge polarization, arises due charge transfer between atoms in response to change in the overlap integr als between neighboring atoms. This exchange charge mode is explicitly modeled in COMB via the charge transfer scheme. No other modes are explicitly included since the COMB model relies on spherically symmetric atom centered density distributions. As an example, the calculated polarizability tensor for molecular oxygen is anisotropic. Polarization parallel to the bond of the molecule is comprised of both short and long range contributions which in COMB are modeled with a redistribution of charge. Polariz ation perpendicular to the bond is only comprised of long range modes. These modes are not modeled in COMB which gives a polarizability of zero in all directions running perpendicular to the bond. Calculated values for the polarizability tensor determined via coupled cluster theory (CCSD(t))87, 88 have a polarizability of zz=2.03 along the bong and xx= yy =0.73 perpendicular to the bond for the molecule. The COMB model misses both the perpendicular response and quite possibly
78 the mitigating effects of the short range modes on the parallel response although the latter e ffects may b e included indirectly in the dynamic charges equilibration. One solution to capture the other polarization modes is to place a point dipole on each atom similar to the fluctuation chargefluctuating point dipole model by Sterne et al.89 With this model, the dipole moment, is calculated directly from the electrostatic field generate by the atomic charges,q iE the neighboring induced dipoles and any external field, extE j N i j j ij ext q i i i iT E E P r E P 1) ( (2 136) N i j ij ij qq ij j q ir r r J q E 04 1 (2 1 37) Pi is the polarizability tensor and Tij is the dipoledipole interaction tensor. For atoms, the polarizability tensor is is otropic and reduces to a scalar value. Induced dipoles calculated in this manner suffer the same instability at close approach as the variable charges.90 Consequently, Tij is employed as a damped functi on that diminishes as atoms overlap. For consistency, the same damping that is used with the Coulombic interactions is applied here: 2 2 2 2 3 02 2 1 1 3 1 4 1ij j ij j r ij ij ij ij ijr r e r r r r Ti (2 139) The additional energy contributions are the dipole self energy, the dipolecharge interact ion and the dipoledipole interactions:
79 i i j j ij i i i q i i i i i i j ij ij i i j j qz ij i i i j j qq ij i i Self i esT E r V J q q J q q V r Q U 2 1 2 2 1 2 12 ) ( ) ( )] ( ), [( (2 138) Since the induced dipoles contribute to the effective field, the point dipole model also requires that induced dipoles are solved self consistently. This leads to an additional iterat ion process during the charge equilibration that can be handled in two ways: either the dipoles are equilibrated self consistently during each charge timestep or the Lagrangian for the system is expanded to include dynamic point dipoles. 9192 The latter approach ensures the time average of the point dipoles, and therefore the Hamiltonian, is conserved, in which case the Lagrangian for the system becomes: ) ( ) ( ), ( 2 1 2 1 2 11 1 2 1 2 1 2 pol N i i i es N i i i N i i Qi N i i iU Q r Q U M Q M r m L (2 139) M is a fictitious mass on the dipoles. The dipole equations of motion can be determined and integrated using the velocity Verlet algorithm as is used throughout this work. In practice, the dipoles converge more rapidly than the charges. A tight convergence tolerance is achieved with one or two iterations of the dipole loop so either method achieves similar efficiencies. In either case, the point dipoles add computational cost and so are only solved for atoms where they contribute a significant portion to the interaction energy, i.e. the interaction between the surface and molecular oxygen. 18.104.22.168 Atom in molecule corrections to the atomic self energies The EE principle has a well defined basis in conceptual DFT,93 a body of work that seeks to explain chemical principles that govern reactivity in terms of DFT. While pursui ng fundamental definitions for such concepts as electronegativity and chemical
80 hardness, that body of work suggests several improvements to consider in EE based potentials. One concept is the variation of atomic hardness when an atom is bound within a mole cule or embedded in an ionic lattice.93 In practice, EE methods are parameterized by either fitting the self energy coefficients to atomic gas phase ionization potentials and electron affinities or by fitting to bulk properties. In the former case, the difference between bulk and atomic properties is compensated by other terms within the potential. Here, the potential is fit to atomic values since the interest is in the interactions betw een atoms and small molecules with surfaces. To capture the change in atomic hardness as an atom is embedded in the bulk oxide, the hardness coefficient is augmented with a correction function that captures the change with its environment: 4 3 2 0, ) 0 ( ) (i ii i i N i j j ij Field ij i i i i i Self iq L q K q q r F J q E q V (2 140) Following the procedure outlined by Troufar et al. ,94 the environmental effect on the atomic hardness is determined by calculating the atomic self energy at the CCSD(t)/aug ccPVTZ95 level in a symmetric field of point charges. The calculation is performed in a field of eight poi nt charges arranged symmetrically around the atom at varying distance. The idea is to determine the effect of a confining potential on the atomic self energy function; the effect is converged with eight point charges. The sum of the point charges compensat es the charge on the central atom so that the net charge of the system remains neutral, thus approximating charge being transferred from the central atom to its neighbors. The Coulombic interaction between the point charges and the central atom as well as between one another is subtracted from the total energy. Thus any change in the atomic self energy function is due solely to the field effect of the point charges on the atomic self energy. Figure 2 2 illustrates the calculated field effect
81 of the lattice on the atomic self energy of oxygen and copper. The effect is found to decay with radial distance as function of r5. A penalty function that captures the change in self energy due to the field strength is fit to the results of the calculation: NN i j ij j J ij j J o j ij Field ir q P r q P q r F5 2 2 3 14 1 ) ( (2 141 ) P1 J and P2 J are adjustable parameters. 2.4.5 Barrier functions The last type of function included in the basic COMB model is barrier functions that restrict the model to regions that are well described by the potential. Functions of this type are not directly related to the chemical bonding but instead limit the potential to conditions where chemical bonding is appropriately described. Ideally, theses functions should not interfere with the calculation but occasionally that is unav oidable. The EE based dynamic charge scheme relies on a continuous self energy function. Such a function can be fit through all charge states until the valence shell limits. For atoms, once the valence shell maximum or minimum is reached the energy as functions of charge has a sharp discontinuity between valence shells. Therefore, the dynamic charge method is only appropriate when charges are limited to the valence shell. In this work, charges are constrained to valence shell with fourth order barrier functions of the form for qi>U iQ: 4) ( ) (i U i Qu i i Barrier iq Q P q V (2 142) Qu is the charge corresponding to an empty valence shell. PQ is a fitted parameter. An identical function is used to penalize partial charges belowL iQ. Charge barriers of this type have also been shown to improve the stability of variable charge potentials.36
82 The second type of barrier function limits the allowed maximum coordination of the atoms. Variable charge schemes tend to favor dense oxide phases with high coordination. An additional barrier function on coordination is needed to destabilize the dense phases such as CuO in the NaCl crystal structure. The coordination correction used in the ReaxFF force field is applied here.16 i Coord i i Coord i i CoordiN N E NV exp( 1 ) ( (2 143) i i iCN CN N (2 144) NN i j ij s ij s ij c iS R r F CN (2 145 ) CN* is the coordination number of the element in the ideal structure, and ECoord and Coord are fitted parameters. Fc is the same as Equation 2135 The correction is only applied for metal oxygen interactions so CNi is the number of M O bonds formed around atom i
83 Figure 2 1. The interaction energy between two charged atoms is plotte d vs. bond length. The interaction between like charges gives positive energies; opposite charges give negative energies. Solid lines are generated using spherical charge densities with =2.0. Dashed lines correspond to point charges. At long separations, the functions are equivalent but begin to deviate near the range of typical bond lengths in solids.
84 Figure 2 2. Self energy functions for O and Cu calculated in a field of eight point charges arranged symmetrically around each atom. Plots correspond to the radial distance ( ) between the point charges and the nuclear center. Calculations are performed at the CCSD(t)/cc pVTZ level of theory.
85 CHAPTER 3 A MULTILEVEL COMPUTATIONAL ANALYSIS OF F LUOROCARBON POLYATOM IC DEPOSITION ON DIAMOND* 3. 1 The C omputational Study of Polyatomic Deposition Reactions The strengths and weaknesses of several computational tools based on various levels of theory are discussed in some detail in the introductory chapter. Such comparisons are best illustrated with a com parative case study that explores the applicability of several levels of theory to a problem of interest. The study presented here compares findings from classical molecular dynamics simulations and density functional theory based simulations of fluorocar bon deposition on diamond. The energies are compared to higher quality multi level post HartreeFock methods which have been shown to yield accurate values for energy in the system of interest. While the study suggests some interesting mechanistic phenom ena about the deposition reaction, perhaps more importantly, it reveals the computational advantages of a classical approach and the consequence of approximat ions and the choices we make while developing classical potentials. Fluorocarbon (FC) plasmas are widely used to chemically modify surfaces in industrial applications. For example, they are used to etch silicon wafers13 and alter the properties of polymeric surfaces in order to increase their chemical inertness and resistance to oxidation.46 The fluo rination process alters the surface free energy and adhesive properties,79 increases the hydrophobicity of the surface,10 results in reduced coefficients of friction,10 and produces higher dielectric constants in electronic devices.11 Similarly, fluorinat ed diamondlike carbon (DLC) films may be synthesized from these Reprinted with permiss ion from Devine, I. Jang, T. Kemper, et al., J. Phys. Chem. C 114, 12535 (2010).
86 plasmas. The resultant DLC films are formed at reduced temperatures compared to those formed from purely hydrocarbon precursors,6 which enables their formation on temperature sensitive substr ates. Thus, the properties of fluorinated surfaces and carbon films can be used for a wide range of applications and, consequently, have garnered significant attention over the past several decades. From an experimental perspective, techniques such as soli d state NMR and x ray photoelectron spectroscopy have established a direct, correlated relationship between products formed during deposition and the final properties of the surface.1214 Furthermore, the products formed on the surface are strongly influenced by the deposition conditions, including pressure, plasma and target temperature, deposition energy, and plasma composition.10,1517 Experimental capabilities have been further enhanced through mass selective ion beam techniques that add a level of dire ct control over the depositing species. With mass selective deposition, the incident particles are separated and accelerated based on the chargeto mass ratio such that the surface products can be reasonably correlated against the mass and energies of the incident ions.2,3,18 What is missing from the above experimental repertoire is a means to explore the dynamic processes and surface chemical reactions that lead to product formation, areas where computational techniques provide increasingly significant i nsight. Often computational studies of surface reactions are constrained to single point calculations within a cluster model. This has been especially true with carbon surfaces, where charge transfer reactions require multi reference methods with sophistic ated treatments for electron correlation, such as complete active space, multi configuration, self -
87 consistent field theory (CASSCF).19 The higher quality results for geometries and relative energies that are achievable with a cluster model have been succes sfully used to explore property relationships, such as the relative stability of different structure reconstructions and relative adsorption energy versus surface terminations.20,21 For example, ab initio methods applied to a cluster model have been used t o map out the mechanisms responsible for chemical vapor deposition (CVD) film growth,22,23 highlighting the critical role of the adsorption processes and hydrogen abstraction from hydrogenated surfaces.24,25 The cluster model has also been used to identify the mechanisms of chemisorption and film growth with radicals on the diamond surface, which often prove intractable by other more efficient methods.2629 Chemically accurate calculations of bond enthalpies are limited to small cluster sizes, which leads t o some deviation from the periodic solid due to a confinement of lattice strain through a smaller volume.29 Findings from mass selected experiments can be readily compared to the results of atomic scale dynamic simulations. As a result of developments i n efficient theoretical models, such as linear scaling density functional theory (LS DFT)30 and reactive empirical potentials,31,32 coupled with dramatic increases in computational power, it is now possible to carry out atomistic simulations of deposition across a range of theoretical rigor. Reactive empirical potentials can efficiently model reactions in systems containing a large numbers of atoms and dynamically simulate a statistically relevant number of deposition trajectories. However, since they are f it to properties of a training set, assumptions are necessarily made when extrapolating to systems beyond this data set. For example, the widely used reactive empirical bond order (REBO)31 potential
88 does not include longrange Coulombic interactions and, s o, is only able to model neutral molecules and radicals. It is fairly well established that during polyatomic ion deposition, and especially during accelerated deposition processes, both cationic and neutral radical species are present and participate in t he reactions. Simulations of these processes using the secondgeneration version of the REBO potential31 must assume that there is no difference in reactivity between cationic and radical states of the same molecule. This assumption may be reasonable at hi gh deposition energies, where reactions predominantly involve fragmentation.16 At lower deposition energies this assumption needs to be rigorously verified and its limits quantified. The middle ground between predictive ab initio methods on small clusters and empirical simulations of large periodic systems is increasingly being filled with DFT, and especially DFT based molecular dynamics (DFTMD). In particular, linear scaling DFT methods may be applied to systems containing several hundred, or even thousand, atoms using modern multi core processors. The SIESTA method and code is a DFT code with linear scaling capabilities, where gains in efficiency are, in part, due to the reduced basis sets .96, 97 This approach may lead to inaccuracies especially when simulating the interaction between a small polyatomic ion and a periodic solid sy stem where basis set superposition error (BSSE) is likely. The numerical basis sets employed in the SIESTA method mitigates BSSE to some extent since the density of the neutral species is well described. Furthermore, the scalability of the SIESTA method requires a basis set that is truncated in spatial extent, which may lead to inaccuracies in the description of surface reactions at extended bond distances.
89 This work explores the surface reactions that occur during polyatomic FC ion deposition on hydrogenterminated diamond surfaces using several different computational methods. In particular, classical MD simulations with the secondgeneration REBO potential and DFTMD simulations are used to explore differences in reactivity between cations and free radi cals in the ion beam. The results are compared to higher level quantum chemical calculations performed on small clusters with the intent of verifying the predictions of the MD simulations. The details of the simulations and each level of calculation used here are described in Section 3. 2. The results for each level of theory are presented in Section 3 .3, followed by a discussion in Section 3 .4 that compares trends observed across all levels of theory. The conclusions of this work are given in Section 3 .5. The simulations in this work were actually performed over several years by a number of coauthors. As a good example of the computational cost of first principles simulations, the DFTD simulations were initially started by Inkook Jang in 2005. At the tim e, a single simulation of one deposition event required nearly six months of processing time to complete using half of the most powerful machines available to us at the time. After several years, Inkook managed to finish the 10 runs in an initial test set Since that time, Donghwa Lee ran the remaining DFTMD simulations presented in the text. Since the initial work machines are much more powerful and the DFTMD algorithms have been substantially optimized so now one trajectory requires only a month on one multi core computer. Inkook also began the classical MD simulations but left the research group before their completion. Travis Kemper, finished that simulation work. All data analysis and the first principles calculations are the work of this author.
90 3.2 Computational Details 3 2 .1 Classical m olecular d ynamics A sampling of the reaction space is conducted using classical MD utilizing the second generation REBO98 potential for carbonbased systems that has been extended to include fluorine.98 Three hundred MD simulations are performed as si ngle, isolated trajectories targeted over a small active region in the center of a larger periodic slab representing the diamond (111) surface. An 11,520 atom supercell is used to model the hydrogen terminated diamond (111) surface. Periodic boundary condi tions (PBCs) are applied in two dimensions within the plane of the surface. The active region contains 128 carbon atoms within an area of 10 2 on the diamond surface. A small active region is used to better enable comparisons to the more computationally i ntensive DFTMD approach discussed in Section 2.2. The deposition energy in the classical MD simulations is 50 eV, which is in a hyperthermal range where smooth FC film growth is observed experimentally4 but where effects due to the excited state potential energy surface are minimized.99 The trajectory is normal to the surface with random variation across the active region. The depositing radicals are also randomly rotated relative to the surface. However, the radicals have little initial internal kinetic energy. The classical simulations are performed at constant volume and temperature (NVT). The system temperature is maintained at 300 K via a Langevin thermostat.100 No frictional force corrections are applied in the active region so that the thermostat does not interfere with the dynamics of the reactions. A Verlet algorithm is used to integrate the equations of motion with a time step of 0.1 fs. Finally, the simulations are carried out for 10 ps, which is sufficiently long for stable products to form.
91 3 .2.2 Density f unctional theory molecular d ynamics Deposition simulations are repeated using DFTMD, which allows us to examine the deposition of both cations and radicals, albeit with a significant increase in computational expens e. Specifically, the DFT implementation within the SIESTA97 methodology and program is used. This approach improves the efficiency of traditional DFT MD simulations for large periodic systems.101 The SIESTA met hod is still subject to same the errors associated with DFT in general, including the approximations present in the exchange and correlation functional, basis set incompleteness and BSSE, over stabilization of partial charges, and static correlation errors to name a few.29 Consequently, an additional objective of this work is to verify the applicability of SIESTA DFT MD for the simulation of FC deposition on diamond. While the SIESTA method has improved efficiency compared to planewave based DFT implementations, it is still computationally expensive. Consequently, only 20 DFTMD trajectories are considered using a supercell composed of the 10x10x13 active region used in the classical MD simulations The supercell contains 128 carbon atoms with hydrogen terminated dangling bonds on the surface (see Figure 1a). In this case, the PBCs are applied in three dimensions, with 15 of vacuum added in the  direction. Two sets of 10 DFTMD simulations ar e carried out; once set under neutral conditions, which corresponds to the deposition of CF3 radicals, and one set with a positive charge applied to the system (that is compensated for using a uniform neutralizing background charge), which corresponds to t he deposition of CF3 +. In each case the incident particles are deposited under the same conditions as during the classical REBO MD simulations (i.e. incident kinetic energies of 50 eV on a trajectory
92 normal to the surface). The simulations are again carried out within the NVT ensemble. However, in this case the temperature is maintained at 300K using a NosHoover thermostat applied to all the atoms.102104 Again, a Verlet algorithm is used to integrate the equations of motion with a time step of 0.1 fs. Simulations are performed for a total of 2 ps, at which point stable products have formed. The nucleus and core electrons of each atom are represented with a Troullier and Martins105 type norm conserving pseudopotential. The valence electron KohnSham states are expanded in a linear combination of atomic orbital type basis sets.96 The basis set parameters for C and F are optimized to the relative energies of CF3 and CF2 radicals and their cations when compared to unrestricted coupled cluster calculations, uCCSD(t)/cc pVTZ.88, 95, 106A double zeta polarized basis (DZ2P) is employed where the polarization functions are explicitly included and determined accordi ng to soft confinement parameters. The optimized basis set is polarized with a d orbital on C and F. The soft confinement potential is set to 100 Ry (1360 eV) with an inner radius of 0.95 times the outer hard cut off. A split norm fraction of 0.25 is used in the doublezeta construction for C and F, while the larger value of 0.5 is employed for H due to its greater variation in effective radius as a function of environment. A mesh cutoff of 250 Ry (3400 eV) is used to determine the fineness of the auxiliary grid basis for the expansion of the density. The above basis set is a compromise between accuracy and efficiency, where a rather small basis set is necessary for computational efficiency in the periodic solid. This leads to a lower quality description of single molecules and, potentially, a significant BSSE when simulating surface reactions However, because the basis functions are
93 generated from the numerical solutions of the atomic problem, the intrinsic BSSE is generally found to be smaller than for Gaussian basis sets, provided the radial confinement is not too great. The PBE functional is used to determine both the exchange and correlation contributions to the total energy in all the SIESTA DFT simulations and spin polarization is used.107 The geometric parameters for CF3, CF3 +, CF2 and CF2 + are presented in Table 1. It can be seen from the table that the DFT method with this basis set overestimates bond lengths by about 2.7% for both CF3 and CF2 cations and radicals, while it underestimates bond angles and improper torsions by 0.7%. The bond lengths and angles are marginally worse for CF2 and CF2 +, with the largest deviations from reference calculations occurring for the CF2 radical. 3 .2.3 Enthalpy c alculations based on a c luster m odel The heat of reaction at 0K ( H0K) is calculated for smaller clusters as a benchmark with which to compare all the methods considered here, and to add some insight into the trends predicted in the classical and DFTMD simulations. Standard references for the H0K are calculated using the multilevel method G3MP2B3.108 The computational cost of this method limits the total system size that can be examined to CF3 interacting with a ten carbon atom cluster (adamantane) that is illustrated in Figure 1b. These calculations are performed using the GAUSSIAN03 ( G03) software package.109 In order to compare these calculations to the findings of the REBO MD and DFTMD simulations, a series of representative reactions are extracted from analy sis of the reactions observed in the simulations. The H0K for each of these representative reactions is then determined using both the G3MP2B3 within G03 and the DFT approach. The DFT calculations performed on the clusters use the same basis set and
94 functional combination as during the DFTMD simulations for both the geometry optimization and energy calculations. The REBO potential is parameterized without corrections for zero point vibrational energy ( ZPE). Therefore, to compare results between methods no corrections for ZPE are made in this calculation. The G3MP2B3 method uses frequencies determined via B3LYP/631G(d) and scaled by a factor of 0.986 to calculate the zero point vibrational corrections.60, 110 This contribution is, therefore, subtracted from the G3MP2B3 enthalpies to provide a direct comparison to REBO and DFT results. Calculations performed with an adamantane substrate do not capture all the energy contributions to the reactions with periodic solids. Ideally, the clus ter should be large enough to encompass the entire bonding region and maintain sp3 hybridization with carbon carbon bonds around the bonding carbon atoms. The ideal cluster size to represent the diamond (111) surface was determined previously by Brown et al.58 In that work, geometries were nearly fully converged within a 22carbon atom cluster (C22), which is illustrated in Figure 1c. This cluster geometry maintains the diamond s tructure around the central carbon atom without freezing the terminal hydrogen atoms, which, consequently, minimizes contributions due to lattice strain to the total enthalpy calculation. This provides a more direct comparison between the various methods s ince the longrange forces vary between DFT and REBO. The multilevel G3MP2B3 method is computationally too expensive to run efficiently with the C22 cluster. Density functional methods, however, are tractable with a decent basis set. Consequently, several hybrid functionals that have shown promise in calculating heats of reactions involving radicals111 are evaluated against the G3MP2B3
95 results for adamantane. In each case, the B3LYP/631G(d) geometries are used. The reaction enthalpies are then calculated using the larger 6311G(2d,p)112 basis set. The specific functionals that are tested along with B3LYP include the BMK113 functional, which has been found to yield promising results with free radical reactions,111, 114 and the B98115 functional, which has been found to produce reasonable bond dissociation energies for a large series of reactions and several molecular properties.114, 116 In all the hybrid energy calculations using DFT with a G aussian type basis sets, the integration grid is set to ultra fine (99,590 pts.) and Gaussian specified tight convergence limits are used in both the self consis tent field (SCF) iteration for energies (106 Hartree) and geometry optimizations. 3. 3 Simulation Results 3 .3.1 Classical MD s imulations To sample the reaction space during radical deposition, 300 independent classical MD trajectories are performed. The si mulations reveal that reactions during deposition are influenced by both the orientation of the deposited radical and the impact site on the diamond surface. Since the CF3 radical has a pyramidal geometry, it reacts differently depending on whether the C atom or F atoms contact the surface initially. Reactions that involve the breaking of a C F bond are only predicted when the F atoms of the FC are oriented towards the surface so that they make initial contact with a terminating H atom. The number of C F bonds broken during initial impact depends on the number F atoms that make contact with H during the initial impact. When the radical approaches the surface with F atoms pointing towards the surface, fragmentation is likely to occur. Likewise, when the radi cal is oriented with the C atom pointing toward the surface, fragmentation reactions are predicted to occur less frequently. Most impacts that
96 produce reactions involve only one F atom contacting the surface and so lead to reactions involving the breaking on one C F bond. The location of the impact site on the surface also influences the reactions that occur. These impact sites can be categorized as an atop site, which is labeled A in Figure 3 2 and corresponds to a direct impact on a surface C atom or its terminal H atom, a subsurface site that corresponds to the recessed surface carbon atom, labeled S in Figure 3 2, and a hollow site, labeled H in Figure 3 2. Each site has equal probability on the surface and, since the trajectories are randomly distribut ed across the surface, trajectories are equally distributed among all three sites. However, because the depositing FC is larger than the reaction sites, a majority of trajectories (60%) involve impact at the A site since it is the most prominent site on the surface. Likewise, the S site is more accessible than the H site, which results in the second highest incidence of impacts occurring at the S site (28%). Impacts that predominantly involve the H site are restricted to fairly specific orientations of the depositing molecule which are predicted to occur only 12% of the time. The impacts on the A sites are reactive and yield new products 80% of the time. In particular, of the A site reactions, 49% abstract a H atom from the surface and 45% result in CF3 bein g bound to the surface with freely dissociated H. The impacts at the S site produce new products 32% of the time, with CF2 and HF being the likely to occur. The CF2 molecule is bound to the surface in about 50% of the reactions in which it is produced. These results are summarized in Table 3 2. The dominant trend revealed by REBO MD simulations is that the A site is significantly more reactive than the other surface sites. Secondly, all the products
97 formed from FC impact on this site involve the removal of H. This is consistent with the accepted mechanisms for growth of diamond films in which radical addition preferentially occurs at the A site following abstraction of atomic H from the surface.13, 117, 118 3 .3.2 DFTMD s imulations Ten trajectories are simulated with the DFTMD approach at both neutral charge and with a positive (+1) charge applied to the system for a total of 20 DFTMD simulations. The results for each charge state are given in Table 3. The orientation dependence that is observed during REBO MD simulations hold true in the DFTMD simulations as well. In particular, reactions that involve the breaking of a C F bond are more likely to occur when the F atom of the deposition species is oriented towards the surface. However, in the DFTMD simulations the FCs appear to be significantly more reactive than in the REBO MD simulations, with reactions occurring during every deposition. In comparing the DFTMD results for the two charge states conside red, we find qualitative similarities and differences in the reactivity of the cations and radicals despite the statistically low number of trajectories. For example, they clearly indicate that the CF3 + cation is able to oxidize the diamond surface. The charge is applied to the system as a whole. As the cation approaches the surface, a Mulliken analysis indicates depletion of electron density from the surface carbon atoms and a neutral charge on the incident FC. The geometry of the FC cation also changes to the pyramidal configuration of the radical species indicating it has been reduced. This suggests that both FC species initially interact with the surface as radicals, which leads to several similarities in the observed reactions. For example, both species tend to form CF2 as a product with
98 similar probabilities (40% for the radical and 30% for the cation). They also both readily form HF when an F atom of the FC directly impacts a surface H atom. When significant fragmentation occurs, both species are able to add to the surface as CFx, where x When the cation oxidizes the diamond slab the surface appears to be less reactive than the neutral surface. This is illustrated by the observation that in 20% of the cationdeposition trajectories the only reaction observed is charge transfer to the cation without any bond breaking, whereas, in the radical deposition simulations, bonds are broken during every trajectory. Additionally, 30% of the cationdeposition trajectories result in H being dislodged from the surface without any occurrence of new bond formation; the pure dissociation of atomic H is observed in only 10% of the radical simulations. In 20 % of the cationdeposition simulations, two H atoms are liberated and combine to form H2. The cations also tend to fragment more extensively producing free atomic F and H species and unbound FCs. In two cationic simulations (20%), atomic F is predicted to migrate to the site of an under coordinated C atom on the diamond surface and form a terminating C F bond. The same reaction is seen just 1 % of the time in REBO simulations and is not observed during DFTMD simulations of radicals. Instead the radicals show a stronger propensity to form bonds with the C atom of the CF3 resulting in either new compounds, such as HCF2, or surface bound FCs. In particular, such reactions are predicted in 50% of the radical deposition simulations versus 20% of the cationdeposition simulations. The statistics of the DFTMD simulations are limited due the small number of trajectories explored and the large number of different reactions predicted. This is a necessary limitation of the method due to its computational cost. However, the results
99 do qualitatively indicate differences between species that can be further explored with calculations using methods of higher fidelity, as was done here. 3 .3.3 Determination of H0K for r eactions on a damantane (C10) The calculations are initially performed on the adamantane cluster illustrated in Figure 1b. Since the reactions of interest involve the A carbon site, calculations are performed on the flat surface of the cluster, as oriented in the figure. The products of radical deposition are illustrated in reactions 16 in Table 4, while the cationdeposition products are summarized in reactions 714 within the same table. Additi onally, reactions 7 and 8 compare different charge transfers that involve the formation of CF2, and reactions 1113 explore different charge transfers involved in H dissociation from the diamond surface. The most energetically favorable radical reaction at the G3MP2B3 level is the attachment of CF2 to the surface with the liberation of HF. The formation of HCF3 is less favorable by 0.34 eV, but is still exothermic. Attachment of CF3 to the surface is exothermic as well, but with an enthalpy that is 0.53 eV higher than the attachment of CF2. The reaction enthalpy for the formation of unbound CF2, reaction 5, is surprisingly high at 2.31 eV considering the prevalence of CF2 in the DFTMD simulations. The 6.12 eV difference in energy between reactions 4 and 5 i s the formation energy of HF. However, these are distinct reactions that are both observed in the DFTMD simulations. In the case of reaction 5, HF is formed during the initial impact. In reaction 4, H and F are dissociated separately by the initial impact Clearly, the hyperthermal deposition energy is driving these unfavorable reactions in the simulations. Similarly, the dislocation of H from the adamantane molecule (Reaction 6) requires 4.64 eV
100 suggesting that this reaction is also driven by the initial kinetic energy of the depositing species. Reactions 714 reveal the most energetically favorable cation reaction to be the formation of the unbound CF2 diradical with a reaction enthalpy of 0.19 eV. The liberation of H from the surface is less favorable by 1.98 eV but also reasonable considering the incident kinetic energy of 50 eV. The large positive enthalpies for reactions 12 and 13 indicate that H departs the surface as atomic H and leaves a cation site on the surface of the diamond cluster. Similarly, the large enthalpy for reaction 7, the formation of CF2 +, shows the formation of the diradical with a cation site in the diamond cluster is more favorable by 5.15 eV. The large positive enthalpies for the bound products in reactions 9 and 10 suggest a dec reased probability for these reactions. For comparison, the calculation of reaction enthalpies is performed with several functionals, as well as with the SIESTA method, and the findings are given in Table 4. The results of calculations using SIESTA with t he same functional and basis set as was used in the DFTMD simulations agree with the results of calculations performed at the G3MP2B3 level with a mean absolute error of 0.36 eV for both radical and cation reactions. The most significant deviations occur for the attachment of the CF3 radical to the surface (reaction 2), where the SIESTA enthalpies are positive and 0.38 eV higher than the G3MP2B3 results, and for the oxidation of the diamond cluster (reaction 14), where the SIESTA DFT enthalpies are negativ e while all other methods give positive values. The B3LYP and B98 results are slightly more accurate with nearly identical
101 absolute errors of 0.23 and 0.22 eV, respectively. The BMK functional has the lowest absolute error with 0.10 eV and so is used exclusively in subsequent calculations. 3 .3.4 Determination of H0K for r eactions on C22 Due to the small size of the adamantane cluster, the results are not directly comparable to those for the periodic diamond slab. This is because of the contribution to the reaction enthalpy due to lattice strain. In addition, the terminal H atoms on the neighboring A sites in the periodic slab are sufficiently close to form hydrogen bonds with a bound FC, such as in reaction 2. These secondary interactions are missing in calculations performed on an adamantane cluster. For better c omparison to the periodic slab, reaction enthalpies for the same reactions are calculated with the C22 cluster illustrated in Figure 1c. In this case, the energies are determined using the BMK functional for exchange and correlation contributions. We find that a large, polarized basis set is necessary to describe the bonding between the CF3 radical and the diamond cluster and to minimize the BSSE. Unfortunately, the cluster size demands that a mixed basis set must be used in the bonding region; the reacting species, surface H atom atoms and the internal C atoms are described with the 6311G (2d,p) basis, while the external C atoms and terminal H atoms are described with the 631G(d) basis. The cluster naturally has a small nonzero dipole moment oriented int o the surface along the surface normal, which the mixed basis set slightly enhances. The effects of the surface dipole are minimized in the enthalpy calculations since it is present in both the reactants and products, but we should expect some distortion of the results, especially when the products involve a charged species on the diamond cluster surface.
102 The calculated H0K values for reactions on the C22 cluster are reported in Table 5. For the radical reactions, once again the formation of HCF3 is favorable. Surprisingly, the binding of the CF3 radical to the surface is exothermic by 0.92 eV with the BMK functional, which is larger by 0.95 eV than the results on adamantane. Reaction enthalpies for fragmentation into CF2 on the C22 cluster differ from the results from the adamantane calculations by only 0.02 eV, indicating that lattice strain due to the formation of a ra dical site in the cluster is nearly relaxed within the adamantane molecule. This finding agrees with earlier studies of radical formation on the diamond (111) surface, which found that the strain effect of radical formation is localized.119 Similarly, trends in the cationic reactions match those of previous calculations; the formation of the CF2 diradical is still exothermic and the binding of the cation to the surface is highly endothermic. A significant deviation occurs for reactions that result in the formation of a cation site in the diamond lattice. The enthalpy for reaction 7 drops by 0.46 eV and for reaction 11 by 0.41 eV due to a more relaxed lattice strain in the larger cluster, and stabili zation of the surface charge by the surface dipole moment. Reaction 14, which considers charge transfer from the diamond cluster to the FC, is more favorable as well, apparently due to the increased stability of the positively charged diamond cluster. Calc ulations are also performed on the C22 cluster using the SIESTA DFT approach and with the REBO potential, the results of which are also reported in Table 5. These results concur with those of the BMK calculations with a mean absolute deviation of 0.05 eV f or REBO and 0.47 eV for SIESTA. Simulations using REBO suggest that the addition of CF3 to the surface is the predominant addition reaction. However, enthalpy
103 calculations predict that the addition of CF2 is energetically more favorable. The most significa nt deviation in the SIESTA based calculations is found for reaction 14 in which the oxidation of the diamond cluster by CF3 + occurs. The SIESTA DFT calculations consistently predict lower reaction enthalpies when a positively charged diamond cluster is for med, which explains the significant decrease in enthalpy for reactions where this occurs. 3.4 Discussion of Results The DFTMD simulations and the enthalpy calculations on both clusters indicate a difference in the reactivity of cations vs. radicals; cations more efficiently remove H from the surface and radicals more readily form compounds. Calculations where the diamond surface is modeled with an adamantane cluster show that formation of a C C bond between the depositing molecule and the diamond surface i s more favorable with the CF3 radical and is energetically unlikely by 4.5 eV with the CF3 + cation. The only energetically favorable reaction, in terms of enthalpy change, that is predicted by simulations of cation deposition is where the positively charged FC oxidizes the diamond target to form the CF2 diradical and a cationic site in the carbon cluster (reaction 6 vs. reaction 7). The reaction appears to be driven by the reduction of the energetic CF3 + cation. The formation of a CF2 + cation as a result of the collision requires 5.18 eV more energy than the formation of the CF2 diradical. The ionization potential of the diamond slab, which is 5.58 eV for the C22 cluster, is offset by the large favorable reduction potential of the CF2 radical. The differenc e in reactivity between charge states is attributable to the cation oxidizing the diamond surface. A consistent trend through all levels of theory applied here is that the formation of molecular products during the initial impact, whether they
104 are surface bound or free, is enhanced with the formation of HF on the surface. With the CF3 radical, HF and CF2 are readily formed whenever an F atom directly impacts an H atom. The oxidized surface inhibits the formation of HF and instead favors the dissociation of unbound H from the surface. HF does form during cation deposition, but more frequently it forms as a secondary product rather than during the initial impact. Similarly, H2 is predicted to be a common secondary product following the recombination of two free H atoms. This overall difference in reactivity, which is initially somewhat surprising given the high impact energy considered, is explained if cations and radicals follow different reaction pathways. The difference in the reactivity between the two char ge states has some interesting implications. Firstly, a beam that contains both charge states at impact will have competing reactions that could alter the composition and structure of the products that form. Secondly, the findings from mass selected deposi tion experiments may not directly apply to more typical production techniques with lower concentrations of cations. Thirdly, classical potentials that do not include charge do not fully capture the chemical complexity of these deposition processes. The ac tual effect of cations on film growth is not easily resolved. 60% of the cation trajectories result in the dissociation of atomic H from the diamond surface to produce free gaseous H atoms. The conventionally accepted mechanisms for growth on diamond films from radical precursors, under low energy conditions, propose a process in which atomic H is abstracted from the surface, leaving behind a radical site.118 Radical precursors may then add directly to the surface at the residual radical s ites. The suspected ratelimiting step in that proposed mechanism is the abstraction of H, at
105 which cationic precursors are predicted to be more effective. Furthermore, the presence of the resulting free atomic H may further enhance subsequent film formati on since film growth rates from radical precursors have been shown to vary directly with the gas phase concentration H.120 The distinction from radical based film growth, however, is that the cationic precursors oxidize the surface leaving behind a cationic site. The role of these cationic s ites in subsequent film growth will have to be investigated further to understand their role during growth. The DFTMD and REBO MD simulations support a similar two step mechanism for film growth with CF3 as a radical precursor, where, the most energetical ly favorable initial reaction is the abstraction of H by the CF3 radical. However, the simulations also suggest an alternate mechanism where the FC radical adds directly to the surface. This onestep addition is predicted with almost equal probability to H abstraction in the REBO simulations. The reaction enthalpies for the singlestep addition (reaction 3) show the reaction to be favorable when coupled to the formation of HF. However in REBO the formation of HF occurs less frequently. To further explore this mechanism and the discrepancy among the methods, the transition state for the single step addition of CF2 to adamantane from a CF3 radial precursor is found using the synchronous transit quasi NewtonRaphson algorithm within the G03 suite (QST2).109 The search and optimization are performed using B3LYP//6 31G(d), as was the case in previous optimizations. The starting point positions of the reactants are varied to replicate several trajectories that produce the products in the MD simulations. A diagram of the optimized transition state is presented in Figure 3. Reactants arrive at the same saddle point regardless of whether the CF3 radical is oriented with its C atom up or down
106 relative to the cluster. The activation enthalpy, determined as the difference in energy between the saddle point geomet ry and the separate reactants, is found to be 7.2 eV. The activation enthalpy is smaller than the 50 eV deposition energy by a substantial amount, so the reaction is able to proceed. The present multilevel analysis provides an opportunity to compare and contrast the computational methods considered, especially in a case such as this where different approximations are made within each method. The first question, when comparing cluster models to periodic systems, is how large must the system be to minimize the effects of lattice strain? In the case of radicals, the radical defect in the diamond structure is contained by its nearest neighbors with little change in enthalpy between the 10 and 22 carbonatom clusters. In contrast, the bound products and the cationic defects require larger surfaces for a realistic description. To further test the size effect, single point calculations are performed using the periodic supercell from the MD simulations, the results of which are presented in Table 6. The calculations show comparable values to those on the C22 cluster, which indicates that the strains due to both a cationic defect in the diamond surface and the addition of CF3 to the surface are contained within the C22 cluster. The predictions from the REBO MD simulations and DFTMD simulations exhibit a few critical differences. Firstly, the classical REBO simulations predict a higher prevalence of nonreactive trajectories; only 58% of the total REBO trajectories result in reactions compared to a 100% reaction probabi lity in the DFTMD radical simulations. Secondly, the products formed most often during the REBO simulations involved CF3, whereas in DFTMD the depositing CF3 radical or cation is never intact at the end of any
107 trajectory. Finally, the formation of HF is predicted to occur much less frequently in the REBO simulations than in the DFTMD simulations. The lack of HF formation in the REBO simulations is initially surprising given that the enthalpy calculations with the same potential favor the formation of HF and CF2. This discrepancy arises due to the short cutoff distance for covalent H F interactions within the REBO potential. In particular, the H F interactions are cut off smoothly between 1.4 and 1.8 This consequently limits the formation of HF during the course of a dynamical simulation. This is illustrated in the transition state shown in Fig 3. The H F bond length in this structure is 1.73 which is near the end of the cutoff region of 1.8 ; in this region the H F interactions are severely damped and when just a 4% strain is placed on the bond, the interactions are terminated completely. In contrast, the C F bond length of 1.78 is just outside the cutoff region for that interaction, which begins at 1.8 and terminates at 2.2 Consequently, wit h REBO, the formation of a CF bond to form CF3 is more likely from the transition state presented in Figure3. This explains the greater occurrence of CF3 remaining intact and the reduced prevalence of HF in the REBO simulations. The short range cutoffs in the REBO formalism are necessary to spatially restrict the bond order contributions. This leads to an improved description of ground state properties, but the detrimental effect of these cutoffs in describing transition states here are apparent. The enthalpy calculations presented here do not include corrections for ZPE as a convenient means to make comparisons between methods. This correction is one of several necessary corrections to first principles calculations needed to compare calculated enthalpies t o experimental values. As a qualitative check, ZPE corrections
108 are calculated for the reactions in Table 5. The last column in Table 5 lists values for H0K determined with the BMK functional and corrected for ZPE. The enthalpy contributions due to ZPE are calculated based on B3LYP//631G(d) frequencies scaled with the G3MP2B3 scaling factor (0.986). The results indicate that the ZPE corrections do not affect the qualitative trends of the calculated enthalpies, but do contribute, on the order of 0.25 eV, to the total enthalpy for the reactions. An important consideration regarding the use of ab initio computational methods based on a LCAO basis sets is BSSE, which arises when reactions involve a small molecule bound to a larger cluster. Out of necessity, D FTMD with the SIESTA method requires a modest basis set for efficiency; the resulting incompleteness of the basis set gives rise to a likely BSSE. To illustrate this, a counterpoise calculation121 is performed for reaction 3 on the C22 cluster using SIESTA DFT with the same basis set as in the previous enthalpy calculations. The H 0K for this reaction is underestimated by 0.18 eV or 17.5% due to BSSE by this method, which is on the order of the reaction enthalpy for the addition of C F2 to the surface. Since the error applies mainly to the reactions that produce a bound product, it is significant and may alter the statistical findings of an MD simulation. In the cluster calculations an attempt is made to minimize this effect by using a larger basis set for all interacting atoms, which is more difficult to do in MD simulations. In practice, the BSSE may lead to an over stabilization of bound products through the course of DFTMD simulations, which must be taken into account when assessin g these predictions. Nevertheless, the qualitative trends in the method are correct.
109 It should be noted that the REBO results are not subject to a BSSE error and may be more reliable where applicable. The enthalpy calculations with REBO clearly indicate th at the potential accurately captures the relative energies of reactant and products involved in the radical reactions. 3. 5 Concluding Remarks The results indicate that the classical REBO potential and the DFT method within the SIESTA program are comparabl e in accuracy to DFT methods using commonly used functionals such as B3LYP when they are applicable. Besides the gain in efficiency, classical methods, such as the REBO potential, are not subject to BSSE and other sources of error specific to quantum mechanical based theories. Unfortunately, the REBO potential is currently limited to charge neutral systems, which, as indicated here, pose a serious limitation to accurately describing experimental cationic deposition. The study also shows that the applicabil ity of REBO is limited since it does not accurately reproduce some of the critical transition states involved in the reactions. The current REBO formalism requires a cutoff between the first and second nearest neighbors to reproduce properties of the equil ibrium ground state. This leads to inaccuracies when describing long range interactions such as bond breaking122 or unusual bonding environments such as amorphous systems.123 The DFT MD method within the SIESTA approach is reasonably efficient for a first principles electronic structure method. The absolute errors in reaction enthalpies are on the order of 0.5 eV, based on the present basis set and functional, which limit its applicability to qualitative analysis under highenergy conditions. The main advantage of DFT MD over classical MD is its transferability, especially to various charge states. Undoubtedly, as improved functionals and wavefunction basedtheories are
110 incorporated, and as gains in computational capacity allows for use of larger basis sets, the quantitative accuracy will improve. The multilevel computational approach considered here suggests that the mechanism for the addition of FC radicals to the H terminated diamond (111) surface does not necessarily depend upon H abstraction. We predict that cationic precursors are less likely to add to the surface unless significant fragmentation occurs. Instead, catio nic fluorocarbons preferentially oxidize the diamond surface leaving behind cationic sites, which is likely to significantly influence subsequent film growth.
111 A B C Figure 31 Target systems are used for calculation of reaction enthalpies. A) The periodic slab composed of a 160 atom supercell project in two directions is used for both c lassical and DFT MD simulations. B) The 10 carbon adamantane molecule allows for calculations up to the G3MP2B3 level on a sp3 hybridized carbon system. C) A 22 carbon atom cluster provides the proper orientation of the surface hydrogen atoms and sp3 hybridization around the bonding site.
112 Figure 32 Impact sites on the diamond(111) surface. The atop site (A) corresponds to the top surface C atom; the subsurf ace site (S) is the recessed surface C atom. The hollow site (H) overlays the void in the surface plane. Figure 33 The geometry of the transition state when a CF3 radical reacts with adamantane to produce a CF2 fragment bound to the cluster surface and free HF. Black spheres are C, dark grey is F and light grey is H. Labels correspond to bond lengths in
113 Table 31. Geometric p arameters for CF23 radicals and c ations. Bond Length () Bond Angle (Deg.) Improp er Torsion (Deg.) a SIESTA b CCSD (T) SIESTA CCSD(T) SIESTA CCSD(T) CF 2 1.344 1.297 103.3 104.9 --CF 2 + 1.249 1.215 123.4 124. 9 --CF 3 1.345 1.314 111.0 111.3 124.0 124.9 CF 3 + 1.259 1.230 120.0 120.0 180.0 180.0 aC omputed at the spinpolarized PBE/DZP ( SIESTA) and bCCSD(T)/cc pVT Z levels of theory Table 32. Classical MD p redicted r eactions b etween CF3 and d iamond (111) Reaction Probability (%) Reaction Site n o reaction 42.0 A,H,S Dia Dia C HCF C CF 3 3 13.7 A H C CF C CFDia Dia 3 3 21.7 A Dia DiaC H CF C CF 3 3 10.0 A F H C CF C CFDia Dia 2 3 5.0 A,S H F C CF C CFDia Dia 2 3 4.0 A,S F C CF C CFDia Dia 2 3 2.0 S H CF FC C CFDia Dia 2 3 1.7 S The Dia subscript indicates the periodic hydrogen terminated diamond target. Table 33 DFT MD predicted r eactions b e tween CF3 and the d iamond (111) surface. CF 3 Radical Reactions Prob.(%) CF 3 + Reactions Prob. (%) no reaction 0. 0 Dia DiaC CF C CF 3 3 20 Dia DiaC H CF C CF 3 3 10 Dia DiaC H CF C CF 3 3 10 Dia DiaCFC F HF C CF 3 30 Dia DiaC H CF C CF 23 3 10 Dia DiaC HF HCF C CF 2 3 10 Dia DiaC F H CF C CF 22 3 10 Dia DiaC HF CF C CF 2 3 20 Dia DiaC F H CF C CF 2 3 10 Dia DiaC F H CF C CF 2 3 20 Dia DiaC HF CF C CF 2 3 10 Dia DiaC F H F C C CF 2 2 3 10 Dia DiaCFC H HF C CF 23 10 Dia DiaC CF HF C CF 23 10 Dia DiaCC H F HF C CF 23 10 The Dia subscript indicates the periodic hydrogen terminated diamond target.
114 Table 34 Reaction e nthalpies (in eV) b etween CF3 and adamantane (C10) Reaction SIESTA *B3LYP *B98 *BMK G3MP2B3 1. 10 3 10 3C HCF C CF 0.29 0.25 0.27 0.26 0.22 2. H C CF C CF 10 3 10 3 0.35 0.29 0.17 0.03 0.03 3. HF C CF C CF 10 2 10 3 0.17 0.12 0.16 0.24 0.56 4. 10 2 10 3: C H F CF C CF 8.82 8.35 8.36 8.41 8.43 5. 10 2 10 3: C HF CF C CF 2.54 2.51 2.52 2.58 2.31 6. 10 3 10 3 C H CF C CF 4.45 4.51 4.50 4.60 4.64 7. 10 2 10 3C HF CF C CF 0.56 0.38 0.31 0.29 0.19 8. 10 2 10 3C HF CF C CF 5.14 4.73 4.82 4.96 4.96 9. H C CF C CF10 3 10 3 4.24 4.94 4.98 4.58 4.52 10. F H C CF C CF 10 2 10 3 10.1 10.4 10.5 10.2 10.1 11. 10 3 10 3C H CF C CF 1.36 1.62 1.66 1.72 1.79 12. 10 3 10 3C H CF C CF 8.34 9.16 9.30 9.21 9.19 13. 10 3 10 3C H CF C CF 4.45 4.51 4.50 4.60 4.64 14. 10 3 10 3C CF C CF 0.45 0.05 0.15 0.28 0.63 Mean Absolute Error (eV) 0.36 0.23 0.22 0.10 -*Energies are calculated with the 6311G(2d,p) basis set. Geometries are optimized with B3LYP/6 31G(d). T =0K
115 Table 35. Reaction enthalpies (in eV) between C F3 and C22 Reaction REBO SIESTA 1 BMK 2 BMK+ ZPE 1. 22 3 22 3C HCF C CF 0.28 0.57 0.26 0.20 2. H C CF C CF 22 3 22 3 0.97 1.03 0.92 0.69 3. HF C CF C CF 22 2 22 3 0.11 0.05 0.28 0.11 4. 22 222 3C H F CF C CF 8.41 8.54 8.43 7.98 5. 22 2 22 3: C HF CF C CF 2.58 2.26 2.60 2.41 6. 22 3 22 3C H CF C CF 4.58 4.17 4.61 4.38 7. 22 2 22 3: C HF CF C CF -1.44 0.77 1.05 8. 22 2 22 3C HF CF C CF -4.86 4.97 4.83 9. H C CF C CF22 3 22 3 -4.92 5.54 5.31 10. F H C CF C CF 22 2 22 3 -10.3 10.7 10.3 11. 22 3 22 3C H CF C CF -0.57 1.25 0.92 14. 22 3 22 3C CF C CF -1.53 0.03 0.30 Mean Absolute Error (eV) 0.05 0.47 --1BMK/6 311G(2d,p) energies with B3LYP/6 31G(d) geometries 2Values are corrected for contributions due to ZPE determined via G3B3MP2 T =0K Table 36. Reaction enthalpies (in eV) b etween CF3 and a p eriodic d iamond slab. Reaction REBO SIESTA 2. H C CF C CF 22 3 22 3 0.97 1.03 2. H C CF C CFDiamond Diamond 3 3 0.96 1.03 5. 22 2 22 3: C HF CF C CF 2.58 2.26 5. Diamond DiamondC HF CF C CF 2 3: 2.54 2.28
116 CHAPTER 4 METHODOLOGY FOR ATOM ISTIC SIMULATIONS OF COPPER OXIDATION USI NG A CHARGE OPTIMIZED M ANYBODY (COMB) POTENTIAL 4.1 Current Theoretical Models of Copper Oxidation Historic ally, copper has been used as a model system for the fundamental study of early stage oxidation and oxide film growth on transition metal surfaces.30, 31, 33, 34, 124126 Notable studies of copper oxidation have been published over the past six decades since copper was included in the original work of Cabrera and Mott.127 Over that time, as the resolution of experimental capabilities has progressed, findings suggest an intriguingly complex oxidation mechanis m that is not yet fully understood. Early studies based on macroscale techniques suggest that oxidation proceeds via the formation of a uniform passivating film as illustrated by the CabreraMott model.34, 124 Subsequent electron microscopy experiments indicate a more complicated mechanism that proceeds through island formation, coalescence, and eventually film formation.32, 125, 128 C urrent in situ electron microscopy experiments, which allow for the direct observation of initial oxidation in controlled atmospheres and are resolved at the atomic scale, have further refined our understanding of the oxidation mechanism. General findings that are consistent throughout a range of studies indicate that oxidation proceeds from chemisorption of atomic oxygen to the formation of an oxygen deficient induction layer that facilitates the growth of three dimensional Cu2O islands.30, 126, 129 The current view is that island growth and coalescence is limited by the surface and interfacial diffusion of oxygen,30 which is in contrast to the passive film models of oxidation such as Cabrera and Mott, where cation diffusion is rate limiting.127 While the current body of experimental work strongly supports a het eroepitaxial growth model limited by oxygen surface diffusion,30, 126, 130 the details of the mechanism,
117 especially during the early stages, are not yet clear and are somewhat contradictory in the literature. For instance, i n situ transmission electron microscopy (TEM) experiments indicated a relationship between island nucleation and step edges.128, 131 However, a more recent study conducted at higher resolution found no significant dependence of nucleation and growth rates on surface defects and step edges, but did find a significant increase in growth r ates at grain boundaries.35, 36 Diffraction and electron microscopy experiments consistently indicate that Cu2O islands form with the  plane of the oxide at the interface regardless of the structure of the metal surface upon which it forms.125, 128 Epitaxial growth oriented along the (110) direction of the metal is preferred.125, 128 While the termination of the oxide is invariant with the metal surface, the shape of the islands does vary among copper surfaces.128 More recently, in situ TEM experiments indicate that island morphology also varies with reaction temperature.36, 130 This leads to the promising possibility that if the mechanism of oxidation can be resolved to the point where the dependence of island morphology on reaction conditions is understood that controlled oxidation could be used as a means to direct the morphology of engineered oxide nanostructures. Resolving the oxidation mechanism requires a theoretical method with the fidelity and capability of efficiently modeling the structures of interest. One choice is an empirical potential that can be incorporated into a molecular dynamics (MD) or kinetic Monte Carlo scheme. First principles methods, such as density functional theory (DFT), have the necessary fidelity but are limited by both the size and the time scale of the problem. Because oxide growth is a diffusion limited process that is influenced by strains at the interface,37, 132 atomic scale simulation of the formation of a single small
118 island involves tens of thousands of atoms, which exceeds the current capacity of traditional DFT methods by at least an order of magnitude. Hence, to enable the study of oxidation on a meaningful length scale, highfidelity analytical inter atomic potentials must be used. Simulations of oxide formation naturally fall under the scope of variable charge potentials that can simultaneously model the oxygen molecule, a metallic solid, and the evolving stoichiometry of the oxide as growth progresses. In choosing a potential form for a metal oxide system, a good first guess would be along the lines of an embedded atom method (EAM)13 or second moment based potential such as Finnis Sinclair21 coupled with an electrostatic model. This exact scheme was pursued by Streitz and Mintmire with their electrostatics+ (ES+) model for aluminum and alumina systems,22 which coupled a variable charge electrostatic scheme to a Finnis Sinclair bond order potential. The ES+ model has subsequently been extended with some refinement to several additional oxide and alloy systems.35, 36 In the case of copper oxide there is significant directionality in the bonding for the most relevant phases. Thus, a potential optimized for covalent interactions, such as the Tersoffs potential,11, 12 may be used. Several recently published potentials link variable charge electrostatics with extended Tersoff type potentials. For example, Yasukawa et al. developed one such potential for the covalent Si SiO2 system.17 This potential form was later refined by Yu et al. in their COMB potential to reproduce the SiO2 phase order.18 The ReaxFF family of force fields employs a similar approach16 and has been adapted to several metal oxide systems, which supports the applicability of this potential form to metal oxides.19, 44, 133 A main concern with using a Tersoff based potential is its
119 applicability to closepacked metals. However, as shown by Brenner, the functional form of Tersoffs bond order potential is algebraically similar to the Finnis Sinc lairs potential with the addition of an angular dependence in the bond order term.52 The applicability of this potential form to model metallic systems comes from the work of Iwasaki et al. who extended Yasukawas potential to several metallic systems15 and, more recently, by Yu et al. who developed a modified form of Yasukawas potential for metallic copper based on Iwasakis parameterizations.18 Yus and Iwasakis parameterizations provide a starting point from whic h to develop an analytical potential for the oxide phases. Copper forms pertinent oxides in two different oxidations states. Oxide islands and the initial oxide layer formed during oxidation of metallic copper are composed of cuprous oxide (Cu2O). The oxide forms in a cubic crystal structure with space group m Pn 3 (Part A of Figure 1).134 The structure can be viewed as a face centered cubic (FCC) lattice of cations, with anions occupying of the tetrahedral sites at positions (1/4,1/4,1/4) and (3/4,3/4,3/4), or as a body center ed cubic (BCC) lattice of anions with cations occupying one half of the (1/4,1/4,1/4) positions. This results in 2fold coordination around each cation with a linear O Cu O bond angle. As oxidation proceeds, the oxide phase converts to the monoclinic cupri c ox ide, CuO, that has space group C2/c (Part B of Figure1).135 The linear O Cu O bond persists in the higher oxide although the cation is now four fold coordinated. 4.2 Computational Methodology The underlying potential used in this study is the COMB potential presented in Chapter 2. While the Yasukawa potential and the subsequent COMB formalism can be
120 parameterized to provide a classical approximation of chemical bonding in the ideal ground state structure, additional bond dependent interactions are required to reproduce the phase order of copper oxide and metallic copper. A similar requirement was found necessary by Yu in the developing the COMB potential for SiO2 and Cu.18 In the current potential, dependence of the bond energy on bond angle and coordination resides solely in bond order function that, in turn, gets smaller with incre asing charge. This tends to be problematic for copper oxides, which exhibit directionality in the bonding in both the Cu(I) and Cu(II) oxides. In both CuO and Cu2O, the O CUO bond angle remains 180o with low coordination around Cu. The bond order terms in bij alone are insufficient to stabilize these structures especially with the charge dependence. Instead, additional charge independent angular and coordination terms are added for each bond type. In Yus parameterization for metallic copper, bij has no angle dependence. Instead, that work uses angular functions based on the first, second and third order Legendre polynomials summed over the nearest neighbors (NN): NN i k i j ijk ijk LP NN i k i j ijk LP NN i k i j ijk LP ijk LP iK K K V, 3 2 1 3 2 2 1 2 1) cos 0 3 cos 0 5 ( 1 ) 0 1 cos 0 3 ( 1 cos 1 ) (cos (4 1) KLP are fit parameters. The angle functions in metallic copper stabilize the FCC lattice relative to HCP at nearest neighbor distances. The terms also allow for a better fit for stacking fault energies than what is typically achievable with empirical models. Equation 4 1 can also be used to stabilize the 1800 O C u O bond angle that is characteristic of
121 the two low energy oxide phases and so is incorporated in this work. The parameters, KLP, are specific to the bond type, i.e. CuCu Cu in the base metal. 4 3 Potential Parameterization The main goal of parameter ization is to faithfully reproduce the surface and mechanical properties of the metal and the relative formation enthalpies of accessible oxide phases while maintaining a reliable degree of transferability to possible nonstoichiometric oxides that may ari se during oxidation. To ensure transferability, parameters are determined as a weighted least squares best fit to properties of multiple phases with variation in coordination around the metal. The fitting process is also restricted somewhat in that the pot ential is an integral part of a family of potentials, which are intended to be compatible with one another. To gain some flexibility in the parameterization, values are defined based on the interaction type. This is a deviation from previous COMB potential s, where interaction parameters were determined from element specific values via mixing rules.18 However, in that previous work, correction functions that we re specific to the interacti on type we re added to improve the potential performance. Since any gain in transferability is compromised with the correction terms, there is no perceived benefit in restricting the parameterization to atomic based values. Furthermore, in practice, the mix ing rules are found to be insufficient to describe critical interactions. For example, the angular coefficients for metals are very weak in the metallic phase. However, the oxide is characterized by a linear O Cu O bond that requires a stronger angular term than what is given by the mixing rules. 4 3 .1 Parameterization of atomic and m etallic c opper The electrostatic self energy terms, VSelf for copper are determined as a least squares best fit to first principles calculations for ionization potential and electron
122 affinities. When VSelf is truncated at the second order the values may be derived directly from the first ionization potential and electron affinity provided that the energy of the neutral atom is taken as the zero point. However, we have found t hat higher order terms are required to model energies beyond the first ionization potential, as is necessary for copper with two relevant oxidation states. The reference states are determined at the coupled cluster level of theory using singlet and doublet and perturbative triplet excitations, CCSD(t).88 A correlation consistent triple zeta basis set (cc pVTZ)95 is found to be sufficient to reproduce the first electron affinity and up to the third ionization potential to an accuracy that is within limits for empirical potentials. All ab initio calculations are performed with the G03 computation suite. The starting point for charge independent parameters are taken from Yu et al.18 which in turn are derived from the Iwasaki potential.15 The cubic phase of Cu2O is basically an inter penetrating network that may be stabilized by covalent CuCu bonding.136, 137 The Cu Cu interaction in the oxide requires a longer cutoff than was used by Yu et al. which in turn necessitates a refit of several parameters. The copper parameters at neutral charge are determined as the weighted least squares best fit to experimental cohesive energy and lattice parameter of the FCC ground state. The training data set also includes un relaxed values for the C11, C12 and C44 elastic constants, (100), (110) and (111) surface energies and the vacancy formation energy, V FE. The energy vs. isometric strain for the FCC ground state as well as BCC and natural HCP structures are also included in the training data set. Unrelaxed values as well as the energy vs. strain curves were calculated with D FT using the Vienna AbInitio Simulation Package ver. 4.6 ( VASP 4.6).138 Energies are determined using the
123 projector augmented plane wave method (PAW). 74 The PW91 generalized gradient approximation (GGA) is used f or the exchange and correlation energies. Kinetic energy cutoff is set at 400 eV. Integration over the Brillouin zone is performed over a 10x10x10 Monkhorst Pack kpoint mesh. The energy vs. isometric strain is calculated over a range of 5% strain and fi t to the Rose equation of state139 which yields an energy vs. strain curve in which the bulk modulus is constant at each point. For fitting purposes, the curves are generated with ionic positions fixed relative to the lattice vectors. The natural equation of state, generated with relaxed ionic positions is used as a final check on the parameterization. Values for cohesive energy and elastic constants differ between calculated and experimental values. To rectify discrepancies between data, the calculated energies and the ground state lattice parameter are scaled to the experimental values using the procedure outlined by Mishin.140 The lattice param eters for higher energy phases are not include in the training set. Only relative energies to the ground state structure are considered. The atomic specific parameters for copper are tabulated in Table 41 and interaction dependent parameters are listed in Table 42. Lastly, the bond dependent angular corrections are tabulated in Table 43. 4.3.3 Parameterization of Molecular Oxygen Oxygen parameters for the self energy function Vself(q) are determined as a best fit to the ionization potential and first and second electron affinity of atomic oxygen calculated at the CCSD(t)/cc pVTZ level of theory. O O interaction parameters are determined as weighted least squares best fit to energy vs. bond length curve for O2, O2 1and O2 2with the intent to capture the bond dissociation energy for the different charge states. Reference data for the energies vs. bond length for the anions are also calculated using coupled cluster theory at the CCSD(t)/cc pVTZ level. Lastly the
124 polarizability tensor for O2, O2 1determ ined at via CCSD/cc pVTZ are also include in the training data set. Atomic parameters for oxygen are presented in Table 41 and O O interaction parameters are listed in Table 42. 4.3.2 Parameterization of copper oxide Energies for the oxide phases are fi t to the enthalpy of formation, Hf, at 0K rather than the cohesive energy. This allows for a comparison between phases of different stoichiometry and charge states. Hf is determined as total energy of the oxide phase minus the energy of the reactants is their reference states: O2(g) and Cu(s). Calculations are performed on ideal structures at 0K for comparison to first principles calculations where appropriate. The copper oxygen interaction parameters and charge dependent copper parameters are fit using a weighted least squares best fi t to Hf of the cuprite phase of Cu2O, CuO in the monoclinic ground state and high pressure NaCl and CsCl phases, and the meta stable paramelaconite phase of Cu4O3. The energy vs. isometric strain and the unrelaxed elastic constants for Cu2O are also included in the fitting data set. Additionally, the training set includes the Hf of a series of phases with varying coordination: the anti fluorite of Cu2O and CuO2 as fluorite, cristobalite and cristobalite. As in the case of the metallic phase, formation enthalpies were calculated using DFT using VASP 4.6. Energies are determined using the same functional, energy cutoff and k point mesh as was used for the metallic phase. The equilibrium lattice parameters for the ground state are scaled to experimental values for cuprous and cupric oxide. Energy vs. isometric strain is determined relative to the scaled lattice constant. DFT typically underestimates the Hf for copper oxides with reported values around 1.24 eV
125 compared to experimental values of 1.75 eV .141, 142143 This is mainly due to the over binding of the O2 reference state in DFT. In this work, the energy of the O2 molecule is fit to a higher level of theory, which reproduces the experimental bond enthalpy. With the im proved reference state energy for O2, the calculated Hf for Cu2O of 1.78eV deviates from experimental values by 1.6%. The difficulties with the O2 reference state in DFT also apply to defect formation enthalpies, Hf def 144 Consequently, values for Hf def for point defects and surface energies that provide a reliable comparison to this potential are not available. Hf def for point and planar defects are used a qualitative assessment of the final potential parameters. Finally, charge dependent param eters are fit with the additional constraint that the electronegativity on each atom is equal with 1.0x106 eV/q in the ideal Cu2O structure. This ensures the ground state is a minimum in both real space and charge space. The full list of interaction parameters are tabulated in Table 42 and 43. 4.4 R esults and Discussion 4.4.1 Properties of Metallic Copper Properties of the metallic phase as predicted by the potential are compared in Table 44 to experimental values as well values calculated from first principles18 with Mishins EAM potential140 and Yus previous C OMB potential18. The elastic properties and surfac e energies are fit to values determined with fixed nuclear positions. The values listed in Table 44 are determined with relaxed nuclear positions. The final parameter set reproduces the experimental elastic moduli and the C11 and C12 elastic
126 constants. H owever, the C44 elastic constant is significantly lower than the calculated values. This is a weakness of the potential that was also observed by Yu et al.18 The cohesive energy of the HCP and BCC simple cubic and diamond cubic phases are included in the training set. The values in Table 44 are based on optimized structures using the final parameter set. The simple cubic and diamond cubic phases are fit with fairly low weight which is reflected in the larger deviations from calculated values. However, the correct qualitative trends are maintained for these low coordination phases. The surface energies are determined as the energy per unit area between the 3D periodic bulk and the 2D periodic bulk with vacuum between the surfaces of interest. The relative energies for the (111), (100) and (110) surfaces agree with the values calculated from first principles18 with a mean error of 3.8%. The surfaces take on a slight charge when relaxed with dynamic charge equilibration that reflects the change in coordination at the surface. This is pictured in Figure 4 2 which shows close packed planes parallel to the (111) surface. The charge alternates between positive and negative with each subsequent plane down from the surface and returns to charge neutrality within seven lattice planes. The graph in Figure 4 2 shows the variation in planar average charge density as determined by the potential. Charge returns to neutral about 10 from the surface. As indicated in the figure, the (111) surface takes on a slight negative charge. The values in parentheses in Table 44 are determined with a charge of zero on each atom and no dynami c charge equilibration. The results are equivalent indicating the surface charge contributes insignificantly to the total surface energy.
127 The formation enthalpies for common point defects and the stacking faults in the (111) plane are calculated as a chec k of the potential. The results are listed in Table 44 and compared to published reference values. The calculations for both planar and point defects are performed on a 36x36x36 3 supercell with periodic boundary conditions applied in three dimensions; t he point defects are optimized at constant volume. The planar defect structures are allowed to relax their volume in the direction normal to the plane of the defect using a steepest descent algorithm. The energies are determined with dynamic charge equilibration. The point defect formation energy is determined relative to the reference state, which is defined as the total energy per atom in the perfect FCC lattice according to the following equation:21 Cu Def FE n nE E ) 1 ( (4 2) Here, n is the number of atoms in the defect structure, Edef is the energy per atom in the defect structure, and Ecu is the energy per atom in the perfect FCC lattice. The predicted formation energies for Cu vacancies and the octahedral interstitial agree with experimental values145 and other methods. A calculation of the formation energi es of the three dumbbell interstitials shows the 100 dumbbell to be the most energetically favorable interstitial, with a defect formation energy of 2.87 eV, in agreement with experimental values.145 Trends in the formation enthalpies for the stacking faults reflect trends predicted by the other reference methods. As shown in Table 44, the potential predicts a bigger formation energy for the intrinsic stacking fault than the previous C OMB parameterization.18 When compa red to the C OMB potential, the predicted formation energy is closer to experimental values but differs from DFT values18 by a greater
128 degree. The predicted extrinsic stacking fault energy is improved; however, the unstable stacking fault energy is underestimated by the current parameter set relative to all other referenced methods. The stacking fault energies are directly influenced by the Legendre polynomial corrections. A careful fi netuning of the parameters in subsequent versions of the potential may improve its performance in applications where improved mechanical properties of the metal are required. Here, the angular corrections for the metal are tuned to optimize properties of both the metal and the oxide phase. As with surfaces, the unstable stacking fault and point defect energies take on a slight charge reflecting the change in coordination. The results of the calculations performed without variable charges are listed in parentheses alongside the charge optimized values in Table 44. The fact that the results are equivalent indicates that the effect of charge is very slight. The differences between this parameterization of the COMB potential and that of and Yu et al. are t herefore attributable to the refitting of the parameters rather than the inclusion of electrostatic interaction. 4.4.2 Properties molecular oxygen The self energy functions for both elements are determined as the least squares best fit to calculated values for electron affinities and ionization potentials. The success of the fit is reflected values listed in Table 45. The potential is able to capture the dissociation behavior of oxygen and its anions as shown in Figures 4 3. COMB predicts the lowest energy state of O2 to be the O2 1 anion with a charge of 1/2 on each atom. This concurs with higher level calculations, which predict the lowest energy charge state to be O2 1. The O2 2 is unstable until it dissociates into separate O1 anions, which also concurs qualitatively with ab initio results up to CCSD(t)/augccpVQZ level of theory. The bond enthalpy for the O2 2
129 predicted by COMB deviates from calculated values more than other ions in the training set. The polarizability of oxygen is fit along wi th the other electrostatic parameters. The starting point for the fit is the atomic polarizability calculated to be 0.716 3 using CCSD(t)/cc pVTZ. The fitted value of 0.36 3 underestimates the atomic values but reproduces the components of the polarizab ility tensor of the O2 molecule that are perpendicular to the bond. The calculated polarizabilities for O2 along with calculated bond lengths and enthalpies are listed in Table 46. 4.4.2 Properties Copper Oxide The lattice parameter, formation enthalpy bulk and shear moduli as well as the C11, C12 and C44 elastic constants for Cu2O were included in the training data set. The accuracy of the fit is reflected in the results listed in Table 47. The potential is fit to elastic properties determined with rigid ionic positions; results in Table 47 are with relaxed positions. The C44 elastic constant exhibits a large drop due to ionic relaxation and consequently exhibits a larger deviation from the reference values. For the oxide phases, defect formation enthalpies Hf def must consider the chemical potential of species added or removed from the system142 according to the following equation: i i i Bulk Def Def Fn E E H0 (4 3) Here, is the chemical potential of species i and ni is the number of atoms of species i added or removed from the perfect structure to form the defect. EDef is the energy of the defect structure and EBu lk is the energy of the perfect structure. i depends on the
130 conditions under which the defect is formed such that i = where is the chemical potential of the reference state.142 For Cu2O to form the relationship 2Cu+O=Cu2O must be true. Furthermore, if we neglect the small dependence on pressure of Cu2O, then 2 Cu+ O= Hf Cu2O. At the copper rich limit, where Cu2O is in equilibrium with metallic Cu, Cu is equal to the chemical potential of metallic Cu Cu(s), which in COMB is taken as the total energy per atom in the perfect FCC lattice. At the oxygen lean limit Cu is zero and O= HF Cu2O. 141, 142, 146 For comparison, values for Hf def are determined at the oxygen lean limit for several point defects; the results are listed in Table 47. The relative energies for the oxygen defects, Vo and Oi agree qualitatively with the results of first principles calculations.141, 142, 146 The formation enthalpy for VCu is higher than the experimental reference values. However, here and in the DFT reference values, the Hf def is determined for the charge neutral defect meaning VCu is formed by removing a neutral Cu atom from the ideal crystal. This is unli kely in real systems where an ion is more likely to be removed. A lower energy charged defect most likely accounts for the discrepancy between the experimental and computational results. Hf def for other defects agree qualitatively with the DFT values cons idering that O varies by 0.51 eV between DFT141 and COMB. Furthermore, COMB allows for calculations on larger structures than are accessible with DFT, where the system is fully relaxed and free of strain,
131 In a manner similar to point defect calculations, the surface energies of the oxide consider the chemical potential of the species added or removed from the perfect crystal to form the surface.147 i i in G A 1 (4 4) Here, A is the sur face area and is the surface free energy. In the COMB potential, the calculations are performed at 0K where contributions due to entropy are zero to compare with reported DFT values. Cu and O at the oxygen lean limit are used for comparison to reported values under similar conditions. 147 The results are given in Table 47 for sever al surfaces. The (111) surface is the only stoichiometric surface listed in Table 47. This surface also gives the smallest deviation between methods. The properties of CuO are fit with light weight, as the Cu2O phase is the main focus of the potential Modeling CuO is complicated by the symmetry breaking JahnTeller distortion that characterizes structures containing d9 transition metals.148 The effect can not be modeled with current DFT implementations.144 Classical potentials have had better suc cess in replicating the effect.148 In the classical models, the Cu2+ cation is modeled as an aspherical ion with a distorted neighbor shell. The COMB formalism has two means of replicating asymmetry in the neighbor shell: the Legendre polynomial angle correction and the point dipole model. In this work, the CuO geometry is not fit explicitly. However, the structure was characterized us ing the final parameter set to see if the COMB formalism offers any advantage in modeling these effects. Results generated with potential are listed in Table 48. Without polarization, the potential predicts a monoclinic structure for CuO with a Hf comparable to experimental
132 values.135 The lattice parameters deviate significantly from the low temperature values.135 However, the total volume of the structure is comparable to DFT results as shown in Figure 4 4. Polarization does not improve the results as shown the last column in Table 48. Table 49 lists Hf for several oxide phases. The experimental Hf for CuO and Cu4O3 are lower in energy than Cu2O149 indicating that Cu2O is only stable under oxygen lean conditi ons. This is illustrated in Figure 4 4, which plots the energy per unit volume for the three low energy phases. In Figure 4 4 the volume and energy vs. isometric strain for the CuO and Cu2O structures compare well with values from DFT. 4.4.3 Behavior of O2 on the Cu surface The behavior of O2 on the Cu (100) surface as predicted by COMB is used as assessment of the potential. The enthalpy of adsorption, HAds, and dissociative adsorption of O2 on the Cu(100) surface is determined by placing O2 at several sites and minimizing the structure via a steepest descent algorithm. HAds is determined using an equation similar to Equation 413 for HDef: i i i s Cu Def Adsn E E H) 100 ( ) ( (4 5) E(100) is the energy of the Cu substrate with a free (100) surface. HAds is determined under oxygen rich conditions, where O=1/2EO 2 (g) and O=0.0.142 The upper row of images in Figure 4 6 shows the initial sites select ed for study. The results are compared with values calculated via DFT in Table 410. When the initial position of the O2 molecule is 2 above the metal surface, the molecule relaxes to a local minimum that is 2.18 above the surface with a HAds of 5. 2 eV. The molecule
133 does not dissociate at this point and only a small charge of 0.092 is transferred to the molecule. When the O2 is placed 1.8 above the surface, the molecule relaxes toward the surface and dissociates at some sites but not at others; the relaxed positions are shown in the lower row of Figure 4 6. The energies in Table 410 indicate that the lowest HAds occurs when an O2 molecule is positioned at a bridge site and relaxes to a state where the O atoms are dissociated and absorbed into neighboring hollow sites. This corresponds to the position labeled Bridge 2 in Figure 4 6. The last row in Table 410 lists values of HAds determined with the point dipole contribution turned off. The overall effect of the fairly small, on the order of 0.1 e V However, the final positions are symmetric, where the point dipoles are very small. The largest effect is seen at the least symmetric final positions. COMB values for Eads are compared to values determined via first principles. The DFT results give the same relative trends and predict dissociation to the hollow sites as being most energetically favorable. However the absolute energy values are substantially lower with COMB than as calculated by DFT. This is partly due to the low binding energy of the O2 molecule in DFT, ( 6.3 eV for DFT vs. 5.2 eV in COMB), which effects all enthalpy calculations involving O2 a s a reference state. Oxygen coverage on the Cu(100) surface and the reconstructions that result have been the subject of numerous theoretical and experimental studies. 25, 39, 150170 Experi mentally, at low oxygen coverage at less than 0.3 monolayers (ML) the surface is characterized by microdomains formed in a c(2x2) arrangement, which can be describe as a O atoms molecules occupying hollow sites. With a ML coverage of 0.3 and above, the oxidized (100) surface is believed reconstruct into a missing row 45 ) 2 2 2 ( R X As
134 a qualitative check quarter half, and full monolayers (ML) of O2 molecules are placed on the Cu surface oriented over the bridge sites. The ML covered structures are annealed for 5 ps at temperatures ranging for 0 to 450 K with a 50 K increment. Up to 0.5 ML coverage the O2 molecules dissociate at the hollow sites at 0 K. The 0.25 ML covered surface is stable up to 500 K. However, oxygen begins to migrate to the subsurface layer at 300 K. The 0.5 ML is stable up to 300 K at which point, O begins to migrate to the first subsurface layer of the Cu substrate and the surface begins to reconstruct. With a full ML of coverage, no dissociation occurs at very low temper atures. The O2 molecules remain intact up to 250 K where the surface begins to reconstructs. Half and full ML covered surfaces reconstruct, however, it can not be determined what surface they are moving towards on the time scale of the simulation. Instead, the 45 ) 2 2 2 ( R X missing row reconstructed was annealed under the same conditions and found to be stable up to 500 K. The reconstructed surface at 300K is shown in Figure 4 7. The fact that COMB predicts this reconstructed surface to at least be a local minimum is encouraging as it means the potential can be used for the study of oxidation on the experimentally observed surface over a wide range of temperatures. 4.4.4 Simulations of the Cu2O(111)||Cu(100) interface The metal oxide interfac e is not well characterized in the literature. This is partly due to the graded nature of the most interfaces, which forms from a reconstructed surface. One means of obtaining an atomically sharp interface is through the electrochemical deposition of Cu2O .171 When deposited on the Cu (100) surface, Cu2O forms a film with initial growth in the (111) direction.171 The resulting Cu2O(111)||Cu(100) interface is atomically sharp and suspected to be semi coherent
135 although this exact structure has not been confirmed experimentally. As a final check of the potential, a model of the interface is relaxed and annealed at 50, 100, 200, 300 and 450 K. The model consists of a 35.3x31.4x36.2 3 supercell containing a 20 thick slab of Cu(s) interfaced with a 15.3 thick slab of Cu2O. The periodic boundary conditions are applied in three dimensions; in other words, no vacuum space exists between the slabs, so the periodic image is one of alternating metal and oxide slabs stacked upon one another. The epitaxial relationship between the metal and oxide is Cu2O(111) ] 2 11 [||Cu(100) which gives a lattice mis match of 3.64% along the epitaxial direction. The slab is annealed under constant temperature and pressure (NPT) for 10ps at each temperature. An MD timestep 0.1 fs was used. The temperature is maintained via a Nose Hoover thermostat.103 The system cell is allowed to relax in three dimensions. During relaxation, the symmetry of the supercell is constrained to remain orthorhombic. The interface is predicted to be stable at all the temperatures considered. The first interfacial layers in both phases rearrange beginning at 100 K suggesting either a lower energy epitaxial relationship or reconstruction exist. Figure 4 8 shows the interface after annealing for 10 ps at 300 K, and the view is along the Cu  or Cu2O ] 2 11 [ direction. During the simulation there is no evidence of a phase t ransformation beyond the first layer or atomic diffusion between phases suggesting the potential is predicting the correct ground state configuration. There is also only a negligible charge transfer between phases, as indicated in Figure 4 9, which shows t he planar average charge density with distance from the interface. The figure indicates that charge transfer is limited to the interfacial region of one Cu layer in the oxide and 12 layers in the metal.
136 4. 5 Concluding Remarks Based on results from the s imulations and calculations, a potential composed of a short range bond order potential coupled with variable charge electrostatics appears to be capable of simulating the ground state bonding environments of Cu as both a pure metal and when oxidized to the Cu2O. The potential also captures the bonding behavior of oxygen as a pure element interacting with Cu(s) and when incorporated in the Cu2O oxide. This enables the large scale simulation of processes such as oxidation of the metal surface or phenomena across the metal oxide interface. Less certain is the ability of this potential formalism to model bonding in CuO, where ligand field effects influence the structure. The potential does predict a monoclinic unit cell for the higher oxide with the correc t volume and formation enthalpy. However, the lattice parameters deviate significantly from experimental results. Since the relative energy for the phases are consistent with experimental values, the potential is useful for studies of the initial stages of oxidation up until the formation of CuO. Future work should seek to improve the potential in this regard. Lastly, the simulations performed in the course of this work are only intended to demonstrate the stability and capabilities of the potential formalism. The statistics are not sufficient to draw any quantitative conclusions about the materials themselves. Instead, a more thorough examination of the materials and processes simulated here are left for future work.
137 A B Figure 4 1. Crystal structure s for the two low energy oxide phases of copper oxide. Large grey spheres represent and small black spheres are Cu. Cuprous oxide (A) forms a cubic lattice with a =4.27. Cupric oxide (B) forms a monoclinic crystal structure with a =4.68 b =3.42 c=5.13 and =99.6o. Figure 4 2 Surface charge in metallic copper. The copper surface exhibits negative charge The charge normal to the  plane is shown. Color corresponds to charge with blue= 0.1 and red0.1 charge units. The plot shows the planar charge density distribution normal to the  surface.
138 Figure 4 3.The bond dissociation energy of O 2 anions: Black marks are fit values. Grey marks are calculated (CCSD(t)/cc pVTZ).
139 Figure 4 4. The energy per unit volume of the three lo w energy phases of copper oxide. Filled symbols represent COMB predicted values. Open symbols represent DFT values.
140 Hollow 1 Bridge 1 Top 1 Hollow 2 Bridge 2 Top 2 Initial Final Figure 45. O2 is placed on the Cu(100) surface in the six orientations pictured and relaxed using the COMB potential. The final geometries are shown in the lower row. The small red spheres represent O and the large orange spheres are Cu. Initial orientations are 1.8 above the metal surface.
141 Figure 46. The relaxed reconstructed Cu(100) surface with 0.5 monolayer coverage of O2. The stable structure has a 45 ) 2 2 2 ( R missing row configuration. Color corresponds to charge with red= 1.3e and blue=1.3e. The larger spheres are Cu. The smaller spheres nestled in between rows Cu are oxygen. Note the alternating positively and negatively charged rows on Cu atoms on the Cu surface which help to stabilize the structure.
142 Figure 47. The Cu2O(111)||Cu(100) interface viewed al ong the Cu direction after annealing at 200 K for 10 ps. Color corresponds to charge with red = 1.3 e and blue =1.3 e. The larger spheres are Cu. The smaller spheres are O. The interfacial region where the lattices are disrupted is confined to the f irst Cu layer in the oxide and one to two layers in the metal.
143 Figure 48. Planar average charge density across Cu2O(111)||Cu(100) interf ace after annealing at 200 K for 10 ps. Values are calculated per atom giving the average charge density per atom in each plane parallel to the interface. The zero point in the graph corresponds to the original interface. The oxide phase in on the positive side of the zero point. There is a slight charge transfer between phases at the interface, which is limi ted to the first 2 planes in each phase.
144 Table 41. Atomic and electrostatic potential parameters Cu O (eVq 1 ) 3.768251 4.700782 J (eVq 2 ) 2.966470 5.064537 K (eVq 3 ) 0.515044 2.756183 L (eVq 4 ) 0.257522 0.992188 ( 1 ) 1.476344 3.012029 (q) 0.293153 0.030819 P ( 3 ) 0.335000 0.323757 P J 1 (eVq 3 r 3 ) 0.470698 0.054039 P J 2 ( eVq 4 r 5 ) 1.086271 1.136518 D u () 0.307561 1.628749 D l () 0.000000 0.244020 Q u (q) 2.000 000 6.000000 Q l (q) 6.000 000 2.000000 CN* 2.20000 0 4.200000 Coord 0.250000 0.249569 E Coord (eV) 0.597603 1.934556 Table 42. Bond dependent potential parameters Cu Cu O O Cu O O Cu A (eV) 712.3527 3523.359 598.2983 7 598.2983 7 B (eV) 102.826 14 204.62585 102.51442 102.51442 ( 1 ) 2.712035 5.516839 3.1293 08 3.129308 ( 1 ) 1.467089 2.527568 1.433963 1.433963 0.231055 2.00000 0.231055 0.231055 1.00000 0 1.00000 1.000000 1.000000 M 1.00000 0 1.00000 1.000000 1.000000 c (rad.) 0.00000 0 43.560000 1.739680 2.043622 d (rad.) 1.00000 0 1.000000 1.000000 1 .000000 h (rad.) 1.00000 0 0.220000 0.297973 0.449820 n B 10.0000 0 10.0000 10.00 000 10.00 000 R s () 3.20000 0 2.2000 3.200000 3.200000 S s () 3.60000 0 2.8000 3.600000 3.600000
145 Table 43. Coefficients for L egendre polynomials Cu Cu Cu Cu O Cu O Cu O 1 LPK (eV) 0.073078 0.069535 0.635646 2 LPK (eV) 0.000000 0.000000 0.635646 3 LPK (eV) 0.019678 0.100428 0.000000 *Parameters are equal to zero for all other bond angles. Table 44. Properties of m etallic copper Exp. DFT g EAM f OMB g COMB a ( 3.6 2 a 3.64 3.62 3.62 3.62 Eo (eV/atom) 3.54 b 3.50 3.54 3.54 3.54 Bulk Mod. (GPa) 139 c 140 140 139 139 C 11 (GPa) 170 c 173 173 170 171 C 12 (GPa) 123 c 123 123 123 123 C 44 (GPa) 75.8 c 80 .1 76.2 49 47.7 (10 6 K 1 ) 16.5 b 14.5 16.9 E o Phase Transitions (eV/atom) HCP 0.006 0.008 0.008 0.008 BCC 0.038 0.046 0.017 0.014 Cubic 0.47 0.43 0.49 0.38 Diamond 1.04 1.08 0.94 0.99 H f Point Defects (eV) V Cu 1.28 d 1.27 1.21 1.17 Cu i oct 2.8 4.2 e 3.06 2.41 4.93 Cu i 100 dumbell 2.8 4.2 2.48 Cu i 110 dumbell 2.8 4.2 6.91 Cu i 111 dumbell 2.8 4.2 8.58 Planar Defects (mJm 2 ) (100) 1780 f 1478 1345 1599 1478 (1478) ** (110) 1780 1609 1475 1646 1519 (1519) (111) 1780 1294 1 239 1295 1218 (1218) ISF 61 g 34.3 36.2 44.4 46.3 (46.3) USF 162 210 161 224 105 (105) Twin 24 g 19.2 18.2 45 23.2 (23.2) *Indicates predicted values. All other value are included in the training set **Values in parentheses are determined with charge equilibration aReference 143 bReference 100 cReference 172 dReference 173 eReference 145 fReference 174 gReference 18
146 Table 45. Ionization potentials (IP) and electron affinities (EA) in eV for Cu and O Table 46. Properties of m olecular o xygen *Reference 143 Cu Ab Initio COMB 1 st EA 1.08 1.08 1 st IP 7.43 7.43 2 nd IP 20.3 18.6 O 1 st EA 1.40 1.40 2 nd EA 6.08 6.08 1st IP 13.5 13.5 Ab Initio/ Experimental COMB Bond Length ( 1.21 1.21 Bond Energy (eV) 5.17* 5.17 Polarizability of O 2 ( 3 ) 11 2.0 2.2 22 0.69 0.65 33 0.69 0.65
147 Table 4 7. Properties of Cu2O (cuprous oxide) Experiment First Principles COMB a ( 4.27 a 4.31 e 4.27 E coh (eV/Cu 2 O) 11.3 b 11.4 e 11.4 H F (eV) 1.73 b 1.24 e 1.7 5 Bulk Mod. (GPa) 112 c 112 1 11 Shear Mod. (GPa) 8.15 8.27 C 11 (GPa) 123 c 123 122 C 12 (GPa) 108 c 107 10 5 C 44 (GPa) 12 c 12.1 57.3 Cu Charge 0.54 Point Defects in Cu 2 O (eV) g V Cu 0.45 d 0.28 1.17 e 2. 49 V Cu_split 0.25 d 0.78 1.24 e 2.20 V O 1.55 e 0.95 O i_tet 1.36 1.47 e 0. 71 O i_Oct 1.69 1.9 e 2.26 Cu i_tet 1.47 e 0.95 Cu i_Oct 1.9 e 0.97 Surface Energies (mJm 2 ) g Cu 2 O(100):Cu 1570 f 1510 Cu 2 O(100):O 1070 f 1030 Cu 2 O(110):Cu 1790 f 2380 Cu 2 O(110):CuO 41 7 f 603 Cu 2 O(111) 785 f 783 aReference 134 bReference 143 cReference 175, 176 dReference 62 eReference 141, 142, 146 fReference 147 gDefect formation and surface energies correspond to oxygen lean conditions. Table 4 8. Properties of m onoclinic CuO Experimental COMB COMB +Polar a ( 4.6 8 a 4. 36 4. 33 b ( 3.42 a 3. 74 3. 77 c ( 5. 13 a 5.04 5.05 (Deg.) 99.6 a 97.1 97.1 H F (eV) 1.6 3 b 1.6 8 1.6 8 aReference 135 bReference 143
148 Table 4 9 Hf (in eV) of copper oxides Structure Experimental DFT/Simulation COMB Cu 2 O Cuprite 1.73 a 1.24 b 1.75 CuO Tenorite 1.63 c 1.11 b 1.6 3 Cu 4 O 3 Paramelaconite 4.72 c 5.74 c 5.11 CuO CsCl 0.53 d 0. 31 CuO NaCl 1.24 d 1.22 CuO 2 Cri stobalite 0.30 d 1.58 CuO 2 Cristobalite 0.34 d 1.31 CuO 2 Fluorite 1.29 d 0.49 Cu 2 O Anti Fluorite 1.14 d 0.44 aReference 143 bReference 142 cR eference 149 DIn house DFT calculations Table 4 10. HAds (in eV) of O2 on Cu(100) surface Hollow 1 Bridge 1 Top 1 Hollow 2 Bridge 2 Top 2 DFT 2.0 2.2 2.2 2.9 4.0 2.4 COMB(Polar) 5.7 6.4 6.7 5.6 7.6 6.2 COMB 5. 6 6.4 6.9 5.6 7.6 6. 2
149 CHAPTER 5 C HARGE O PTIMIZED POTENTIALS FOR ALUMI NUM AND ALUMINUM OXIDE SYSTEMS: SUCCESSES, CHALLENGES AND OPPORTUNITIES 5.1 Background and Motivation Aluminum is one the most widely used metals due mainly to its abundance and excellent strength to weight ratio, but also due, in part to the corrosion resi stance afforded by a readily forming protective oxide layer.177 Properties of the oxide layer have long been exploited to enhance the performance of the metal in numerous industrially significant applications such as thermal barrier coatings178, corrosion and wear resistant coatings177, as well as technologically important applications such as microelectronics and catalysis179, 180. In applications that rely on the protective qualities of the coating or where the metal interfaces with another material, performance and failure are often dependent upon interfacial properties and processes. Consequently over the past several decades a s ignificant amount of fundamental and applied research has been devoted to exploring the properties of the metal oxide interface.36, 40, 181188 Given the significance of this interface, it is a reasonabl e undertaking to extend the COMB formalism to model these materials. Crafting a potential for composite interfaces with the efficiency and scalability to emulate lifesize nanoscale devices is not a trivial task. The challenge is further complicated for AlAl2O3 systems since the complex crystal structures of aluminum oxides are not amenable to empirical models.39 While the COMB formalism and similar models are sufficiently flexible to produce a working model of the some aluminas ,40, 189, 190 getting the correct phase order in the ground state requires a well designed potential and a delicate fine tuning of parameters. In the end, all current potentials for this material compromise the accuracy in some properties in
150 favor of others. The final model developed here is no different from existing potentials in that regard. This project began as the initial test project for the development of COMB. At the start, the actual form of the potential was not established, so as a first gues s, three existing potentials were incorporated in an existing in house code. This was a collaborative effort with three researchers: Alan McGaughey coded up the variable charge potential devel oped by Streitz Mintmire, Jiang u o Yu coded the Yasukawa potenti al upon which COMB is based, and Inkook Jang coded a Modified EAM potential. None of the potentials worked very well initially and required significant trial and error to find a reliable potential form. In time, the current COMB model developed as a combination of Streitz and Mintmires electrostatic functions with Yasukawas bond order formalism. This current version would not have been settled upon without going through the effort of coding and evaluating all three potentials. In addition to the COMB code, these three researchers and I have written all the fitting and analysis codes used to develop COMB models for new materials. This is a prohibitively burdensome task, which could not have been completed wit hout the efforts of Alan, Jiangu o and Inkook 5.2 Challenges with Current Theoretical Models for Alumina Experimentally, an oxide layer forms on the metallic Al surface as an amorphous layer of roughly 36 nm thickness.127, 191 The naturally occurring interface formed during thermal oxidation is difficult to characterize at the atomic scale, both computationally and experimentally, since the amorphous interface is graded with mixed stoichiometry over an interfacial region188. Consequently, the majority of work at the atomic scale involves ideal interfaces that represent the key aspects of the interfac e. For example, the most notable experimental characterizations of the solid/solid interface were
151 performed on the interface formed by the deposition of metallic Al on a crystalline Al2O3 substrate rather than the thermally oxidized metal surface.181 This approach was used by Medlin et al. who used high resolution transition electron microscopy (HR TEM) to deter mine the structure and preferred orientation of Al films grown on the (0001) surface of Al2O3.188 The work revealed an atomically sharp interface with three predominant heteroepitaxial orientations, with the preferred relationship of 3 2) 0001 ]( 0 1 10 [ ) 111 ]( 10 1 [O Al Al.This orientation relationship aligns the closed packed planes of the metal with t he closed packed oxygen sublattice in the oxide. The work of Medlin et al. was performed at 200 C; Vermeersch et al. report results from a similar experiment at room temperature where they found that growth of the metal film proceeds through the formation of a suboxide layer, to island formation and eventually, to coalescence and film formation.192 At higher temperature Vermeersch et al observed film growth with the heteroepitaxial relationship of 3 2) 0001 ]( 01 1 2 [ ) 111 ]( 1 1 2 [O Al Al which is the secondary relationships found by Medlin.192 Computational studies based on first principles methods have explored the structure and adhesion of the ideal atomically sharp interfaces. Batyrev et al used DFT to examine the structure and work of adhesion of the close packed interface in the preferred orientation found by Medlin.86 Zhang and Smith looked at variation in the stoichiometry and adhesion of t his same interface with oxygen partial pressure and found that the ideal termination of the oxide depends on the chemical potential of oxygen.185 In a similar DFT experiment, Siegle et al. compared works of adhesion (WSep) in the preferred Medlin interface with several terminations of the Al2O3 (0001) surface and found that when the oxide is terminated with a singl e layer of Al atoms, the
152 calculated WSep corresponds to experimentally measured values.184 More recently Hinnemann and Carter looked at the initial stages of f ilm formation by studying the adsorption of Al and O on the Al2O3 (0001) surface using DFT where they found Al prefers to adsorb at the three fold hollow sites with significant charge transfer from the metal atom to the terminal O atoms.182 The main limitation with first princi ples methods is the computational cost. Wilson and Finnis have shown that prediction of the correct structure of the alpha phase requires a high energy cutoff, at least 600 eV.41 An interfacial study requires a large supercell with sufficient depth in both phases to converge the energy to a consistent point. The c omputational burden restricts the applicability of first principles methods to ideal, coherent crystalline interfaces wi th the preferred orientation of 3 2) 0001 ]( 0 1 10 [ ) 111 ]( 10 1 [O Al Al. First principles based analysis of semi coherent interfaces or other heteroepitaxial arrangements is prohibitively expensive. The development of empirical variable charge potentials that have the computational efficiency to model much larger systems have been pursued over that past decade in response to the limitations in first pri nciples methods. The main difficulty in developing an empirical potential for the Al/Al2O3 interface is finding a suitable potential for the oxide phase. The phase order in alumina has been a challenge for empirical potentials for several decades.39 Historically, aluminum oxides have been modeled as ionic materials with polarizable ionic potentials, such as core shell or breathing core shell models80, 193. Gale et al have shown that these simple polarizable models are unable to predict alumina (corundum) as the low energy phase of Al2O3.39
1 53 Instead, the low energy phase is predicted by these potentials to have the more symmetric Bixbyite crystal structure .39 The two crystal stru ctures are illustrated in Figure 5 1. Part A of the figure shows the hexagonal unit cell of phase with space group c R 3 .134 The crystal is composed of a partially occupied layer of Al stacked in the (0001) direction and sandwiched between two close packed oxygen sublattices. The distinguishing feature of the crystal is short separation distance of 2.66 between one pair of cations associated with each Al atom. The smallest cationcation separation in the oxide is actually smaller than the m etallic bond distance in the FCC metal. Part B of Figure 5 1 shows the cubic Bixbyite crystal. The entire structure, with space group 3 Ia is built up as a 2x2x2 arrangement of subunits that are similar to a fluorite structure consisti ng of a close packed FCC Al lattice with O occupying of the tetrahedral sites.134 The symmetry is dictated by the relative arrangement of the unfilled tetrahedral sites. The local bonding around each atom is more symmetric in Bixbyite than corundum. Corundum has two Al O bond lengths of 1. 86 and 1.97 while the ideal Bixbyite structure has one. The bond angles in the Bixbyite structure are also more symmetric. The Al O Al bond angles are closer to the ideal tetrahedral angle of 109.47o as shown in Figure 5 2 which compares the distribution of the Al O Al bond angles in each material. Bixbyite also has 180o O AlO bond angles, which do not occur in the corundum structure. Gale surmised that the corundum structure may be stabilized with a polarizable potential by including quadrapole m oments in addition to the breathing mode on the anion.39 Around the same tim e, using DFT, Wilson and Finnis found the Bixbyite structure is favored for alumina up to a cutoff energy of 600 eV beyond which the
154 corundum structure is predicted.41 A plot of the difference in the electron density distribution between these two cutoffs reveals that the quadrapole moment of the anions is the distinguishable feature that differs between energy cutoffs.41 The theory was further confirmed by Wilson who compared a series of classical polarizable potentials that sequentially added a breathing mode, dipole moments and quadrapole moments.41 The experiment confirmed that the inclusion of the quadrapole moments of the anion can be used to stabilize the corundum structure. While the quadrapole models successfully predict the correct phase order and appear to be the most transferrable potentials for aluminum oxide, they are limited by a relia nce on fixed, formal charges in current versions.38, 41 This excludes their use from applicat ions containing mixed oxidation states such metal oxide interfaces. Variable charge schemes reproduce some polarization effects, such as exchange charge polarization between bonded neighbors,80 which raises the possibility that a potential of this type may be applicable to alumina. However, variable charge potentials for this system have had mixed success in this regard. The first reasonable model was developed Streitz and Mintmire who coupled a variable charge electrostatic potential with a charge independent Finnis Sinclair potential21 to develop the Electrostatics Plus (ES+) model for Al and Al2O3.22 This model fails to predict corundum as the low energy structure.40, 86 In interface applications it predicts disordered interfaces characterized by rapid inter diffusion of O and Al between phases,22, 187 processes that are not observed experimentally or via first principles calculations.184, 188 Zhou improved the stability of the ES+ model with barrier functions on the charges and by replacing the short range potential with an EAM model fit to various oxygen coordination
155 environments.36 Most recently, Zhang et al. published an extended Tersoff based variable charge potential (ReaxFF) for aluminum alumina systems.40 In that work, the authors found that neither bond bending nor torsion terms are required to predict the correct phase order of corundum. Instead, they found that cor rect structure is predicted when dispersion interactions and a penalty function for over coordination are included in the potential. Zhangs ReaxFF model does not include bond angle dependent contributions and so should describe the same physical interaction as the ES+ model. The difference between these potentials is the inclusion of dispersion interactions and an additional coordination function, which seems to be the key interactions needed to describe the structure of the phase. Wilsons polarizable model also includes a large dispersion interaction between oxygen pairs. The quadrapole moments are a small contribution to the total energy, which appear to introduce a short range angular dependence to the energy. Considering the successes of Zhangs and Wilsons models, the COMB formalism appear to be sufficiently comprehensive to model both alumina and its interface with metallic Al. 5.2 Parameterization of the COMB Potential The COMB formalism used in this work is identical to what was described in Chapter 2. The angular term developed in Chapter 4 (Equation 41) applies here as well. As with the Cu/CuO potential described in Chapter 4, parameters are determined as a weighted least squares best fit to properties of multiple phases with variation in coordination around the metal to achieve a degree of transferability suitable for the simulations of interfaces.
156 The same oxygen parameters that were developed in Chapter 4 are used in this study, as well. However, the O O cutoff distance is extended f rom 2.2 to 3.0 so O O interactions contribute to the energy of the oxide. The dispersion interaction betweens O O pairs are als o included using the following equation since Zhang and Wilsons methods both indicate the importance of these contributions: 6 2 1) (ij VdW j VdW i VdW ijr P P V (5 1) PVdW are fitted parameters specific to each element. In the final parameter set, the best fit is found with zero values for the dispersion coefficients. 5 .4.1 Parameterization of atomic and m etallic Al The electrostatic sel f energy terms, VSelf for atomic Al are determined as a least squares best fit to the 1st, 2nd and 3rd ionization potentials and electron affinity calculated at the coupled cluster level of theory using singlet and doublet, and perturbative triplet excit ations, CCSD(t).88 A correlation consistent triple zeta basis set (cc pVTZ) is found to be sufficient to reproduce the first electron affinity and up to the third ionization potential with sufficient accuracy.95 All ab init io calculations are performed with the G 03 computational suite. The starting point for charge independent Al Al bond parameters are taken from the Iwasaki potential.15 The parameters are fit using a weighted least squares best fit to experimental cohesive energy and lattice parameter of the FCC ground state. The training data set also includes un relaxed values for the C11, C12 and C44 elastic constants, (100), (110) and (111) surface energies and the vacancy formation energy. The energy vs. isometric strain for the FCC ground state as well as BCC and natural
157 HCP structures were also included in the training data set. Unrelaxed values as well as the energy vs. strain curves were calculated with DFT using VASP 4.6. Energies are determined using the PAW method. The PBE functional is used for the exchange and correlation energies. Kinetic energy cutoff is set at 7 00 eV. Integration over the Brillouin zone is performed over a 10x10x10 Monkhorst Pack kpoint mesh. The energy vs. isometric strain is calculated over a range of 1% strain and fit to the Rose equation of state139 which yields an energy vs. strain curve in which the bulk modulus is constant at each point. For fitt ing purposes, the curves are generated with ionic positions fixed relative to the lattice vectors. The natural equation of state, generated with relaxed ionic positions is used as a final check on the parameterization. Values for cohesive energy and elasti c constants differ between calculated and experimental values. To rectify discrepancies between data, the calculated energies and the ground state lattice parameter are scaled to the experimental values using the procedure outlined by Mishin140. The latt ice parameters for higher energy phases are not include in the training set. Only relative energies to the ground state structure are considered. The atomic specific parameters for Al are tabulated in Table 51 and interaction dependent parameters are li sted in Table 52. The bond dependent angular corrections are fit to the intrinsic, unstable and twin stacking fault energies (SF, US andT) as well as, the energy difference between the FCC and HCP phases and the elastic constants Results are liste d in Table 53. Dispersion is not explicitly included in the Al Al interactions since the current potential form cannot model the change in these interactions with charge. Explicit inclusion of dispersion between cation pairs would require the simultaneous fit of both
158 metal and oxide properties, which is considerably more expensive and requires a compromised solution for at least one phase. Instead, dispersion interactions in the charge neutral ground state are assumed to be included in the empirical parameterization on the Al Al interactions. I n the oxide, dispersion interactions are assumed to be smaller between the less polarizable cations, an effect that is captured by variation in the short range att raction terms with charge. 5.4.2 Parameterization of aluminum oxide The Al O interaction parameters, charge dependent Al parameters and bond angles dependent terms are fit using a weighted least squares best fit to the Hf for Al2O3, Bixbyite and Al2O3, as well as, AlO as NaCl. O2(g) and Al(s) are used as the reference states. The unrelaxed elastic constants, bulk and shear moduli for Al2O3 as well as the energy vs. isometric strain are also included in the fit. As in the case of the metallic phase, relative energies of the reference systems are calculated using DFT. Again, the PBE functional is used for the exchange and correlation energies. Kinetic energy cutoff is set at 7 00 eV. Integration over the Brillouin z one is performed over a 6x6x6 Monkhorst Pack kpoint mesh. The equilibrium lattice parameters are scaled to experimental values for Al2O3. Energy vs. isometric strain is determined relative to the scaled lattice constant. The full list of interaction parameters are tabulated in Table 52. Hf def for point and planar defects are used as qualitative assessment of the final potential parameters rather than included in the training data set since these values subject to the choice and accuracy of the reference states.194 Hf def is determined as a small energy difference between two bulk phases. In practice, large variations in these
159 values arise with small changes in the bulk energy, which make fitting to them particularly difficult. All charge dependent parameters are fi t with the constraints that the electronegativity of each site is equal within 1.0x106 eV/q, where q is a formal charge unit, in the phase at 0 K. Furthermore, the electronegativity in the ideal structure is fit to be equal to zero meaning atomic charges correspond to a minimum in the energy with respect to charge without constraints. This is an arbitrary condition where one could expect the actual oxide phase to be either electropositive or electronegative. However, by ensuring that the ideal structure corresponds to a minimum in the total energy, aberrant phases, such as Bixbyite, naturally form at higher energies. This condition is found to be necessary to raise the energy of the Bixbyite phase relative to the alpha phase. 5.5 R esults and Discussi on 5.5.1 Properties of metallic aluminum Properties of the metallic phase as predicted by the potential are compared to experimental and calculated reference values in Table 54. The elastic properties and surface energies are fit to values determined with fixed nuclear positions. Values listed in Table 54 are determined with relaxed nuclear positions. The cohesive energy of the HCP and BCC simple cubic and diamond cubic phases are included in the training set. The values in Table 54 are based on optimi zed structures using the final parameter set. The simple cubic and diamond cubic phases are fit with fairly low weight which is reflected in the larger deviations from calculated values. Qualitative trends are maintained with these low coordination phases
160 The surface energies are determined as the energy per unit area between the 3D periodic bulk and the 2D periodic bulk with vacuum between the surfaces of interest. Relative energies for the (111), (100) and (110) surface agree with values calculated from first principles40 with a mean error of 1.0%. As with metallic Cu, in Chapter 4, surfaces take on a slight negative charge when relaxed with dynamic charge equilibration reflecting the under coordination at the surface. The c harge alternates between positive and negative with each subsequent plane down from the surface and returns to charge neutrality within 10 lattice planes. As with Cu, the surface charge does not make a detectable contribution to the total surface energy. The formation enthalpies for a common point defects and the stacking faults in the  plane are calculated with relaxed ionic positions as a check Results are listed in Table 54 with comparison to reference values. Calculations for both planar and poi nt defects are performed on a 36x36x36 3 super cell with periodic boundary conditions applied in three dimensions. Point defects are optimized at constant volume. Planar defect structures are allowed to relax their volume in the direction normal to the pl ane of the defect using a steepest descent algorithm. Energies are determined with dynamic charge equilibration. The point defect formation energy is determined relative to the reference state defined as the total energy per atom in the perfect FCC latti ce according Equation 41. The predicted formation energies for Al vacancies and the octahedral interstitial agree with experimental values and other methods. 40, 140 A calculation of the formation energies of the three dumbbell interstitials shows the 100 dumbbell to be the most energetically favorable interstitial, with a defect formation energy of 2.50 eV.
161 Trends in the formation enthalpies for the stacking faults reflect trends predicted by other methods as shown in Table 54, with the exception that the potential predicts a large formation energy for the unstable stacking fault. The stacking fault energies are directly influenced by the Legendre polynomial corrections. In this fit, the angle functions also affect the phase order, elastic constants and surfa ce energies The best solution results in a larger value for US. Future refitting may address this weakness and improve the performance of the potential in regards to mechanical properties. 5.5.2 Properties of aluminum oxide The self energy functions in COMB are fit directly to the atomic ionization energies. The accuracy of the fit are reflected in the calculated values listed in Table 55. In previously published potentials, there are two approaches to fitting self energy functions. The method fol lowed by Rappe and Goddard 14and in the subsequent ReaxFF potential40 fits the self energy to second order to atomic ionization values. Other potentials such as the ES+ model and the first generation COMB potential18 fit the self energy term to properties of the bulk phases. Here, the approach of Rappe and Goddard is followed since one intended use of the potential is the study of oxidation, where i nteractions of single atoms with the surface are important. However, the function is terminated at fourth order rather at second order as in Rappe and Goddard. The bulk and elastic constant s are include in the training data set and so are fit directly. T he final values are tabulated in Table 56 in comparison to experimental and first principles values as well as published values determined with ReaxFF and the ES+ potentials. The values for the elastic constants and moduli in Table 56 are determined
162 wit h relaxed ionic positions. Overall, the bulk and elastic constants are comparable to other potential and the reference data. Calculated values for Hf def for charge neutral point defects are listed in Table 56 and compared to published calculations. T he reference values are from Hine et al .194 and are determined using Equation 43 under oxygen lean conditions at T=2371 K and pO2=0.2 atm. Reference states are Al(s) and O2(g). The values in Table 56 are not corrected for o as is the case in Hine et al.194 Zhang et al have shown through first principles calculations that under oxygen lean conditions the preferred termination of the (0001) surface of Al2O3 the interface is with a single Al lay er.185, 186 This observation agrees with experimental data obtained on the clean surfaces under ultrahigh vacuum.181, 195 In the presence of hydrogen or at higher oxygen partial pressures, O terminated surfaces are favored.181 The single Al terminated surface is stoichiometric and so provides a direct comparison between calculations without the need to consider the chemical potential of species added or removed to for m the surface. For these reasons, the (0001) surface energy, (0001), of the single Al terminated surface is used as a check of the potential. The calculated value is 1.82 J.m2. With first principles, functionals based on the local density approximation (LDA) tend to give higher values around 2.2 and GGA functionals give lower values around 1.6. 184, 185 The COMB result falls between these two levels of theory. The COMB formalism stabilizes the corundum structure relative to Bixbyite by 0.74 eV/Al2O3, which is inline with first principles calculations.19 The potential is fit to both the Bixbyite and phases similarly to Wilsons and Zhangs potentials.19, 41 These
163 other phases lie close in energy to the phase. The potential is unable to stabilize the phase relative to Bixbyite as shown in Figure 5 3 which plots the energy vs. volume for several oxide phases. In Figure 5 3, the phase remains continuously above the curve indicating that, with this parameterization, the phase would never form. This is a weakness of the current model, which indicates limited transferability of the potential to other phases. The rocksalt structure of AlO is included in the fit. A coordination correction function (Equation 2129) is required to destabilize this phase relative to the This structure is predicted to be the lowest energy configuration in terms of energy per atom with the ES+ model. Calculations using the ES+ potential withi n our in house code predict this phase to be 0.71 eV/atom lower in energy than the phase, which leads to rapid diffusion of Al atoms from the bulk metal into the oxide during interfacial simulations. A similar effect occurs within COMB if the coordinati on based correction functions are set to zero. Polarizability is introduced via the point dipole model described in Chapter 2. As a check of the effect of polarization, the Bixbyite phase was optimized with polarization both enabled and disabled. When the point dipole model is active, the Bixbyite structure is further destabilized by 0.23 eV/Al2O3, indicating that the point dipoles have some effect on the phase order. The time per iteration is roughly doubled with the point dipoles enabled for the entire system. Since polarization is not necessarily required to stabilize corundum, polarization is disabled in all subsequent calculations.
164 5.5.3 Properties of the Al(111)|| Al2O3(0001) interface To test the reliability of the potential in interfacial appli cations, WSep was determined for the several variations of the 3 2) 0001 ]( 0 1 10 [ ) 111 ]( 10 1 [O Al Alinterface. WSep is a measure of the interfacial energy defined as: Tot Slab Al Tot O Al Tot Slab O Al SepE E E W3 3 2 3 2 (5 2 ) The calculation is performed on a supercell consisting of a peri odic slab of Al metal interfaced with a slab of Al2O3. Periodic boundary conditions are applied in three dimensions so the structure actually consists of alternating films of Al and Al2O3 sandwiched together. Each slab is 40 thick. Previous theoretical studies on this interface conclude that the most stable interface under oxygen lean condition is formed with a single al atom terminating the oxide layer.40, 184186 The same interface termination is used here for comparison to these studies There are six possible variations of the single Al terminated interf ace. The first layer of the Al slab (Al1), can be oriented over three unique sites on the Al2O3(0001) surface relative to the oxygen sublattice. These sites are illustrated in Figure 5 4 where the al slab can be aligned with an FCC (F) or HCP (H) orientation or directly on top (O) of the terminal oxygen layer. Furthermore, as shown in Figur e 5 5 with each configuration of the first Al layer, the second Al layer in the slab (Al2) can be oriented over either of the two sites not occupied by Al1. The combination results in six variations in interfacial structure. The six interfaces were annealed at 50 K for 10 ps and optimized using a steepest decent algorithm. The WSep and the relaxed interfacial separation distance, Do, are compared to other published values in Table 57. The resulting values indicate
165 that the potential predicts very litt le reconstruction of the interface in comparison to DFT.184 The values for both Do and WSep are comparable to unrelaxed interfaces. Additionally, the highes t and lowest energy interfaces are plotted in Figure 5 6. The two images are nearly identical, showing very little reconstruction and some charge transfer between phases. The average distribution of atomic planes normal to the interface for the low energy interface is plotted in Figure 5 7. The distribution indicates very little alteration of the structure at the interface. Average charge density per unit area is also plotted in Figure 5 7. The density plot indicates a slight charge transfer between phases at the interface. The current param eterization predicts less reactivity at the interface than is indicated by first principles calculations Annealing at high temperature does not initiate reconstruction suggesting that the interfacial geometries predicted by first principles based methods lie high in energy than what is predicted by COMB. 5.5 Concluding Remarks This work has demonstrated that the COMB formalism has sufficient flexibility to model aluminum oxides. The main challenge in modeling this material is predicting the correct ground state structure. COMB is able to be parameterized to favor the correct ground state configuration. In the version presented here, less emphasis was placed on the other naturally occurring low energy phase such as The lack of emphasis, in turn, leads to the phase being too high in energy to the point where it does not form. In a similar way, the parameterization may have raised the energy of the DFT predicted relaxed interfaces thus inhibiting any reconstruct ion at the interface. However, further investigation is required before any conclusion can be made in that regard.
166 The potential reproduces the properties of the metal and its surfaces. It also makes use of the O2 parameterization developed in Chapter 4 for copper oxides. Together, both parameterizations extend the applicability of the current potential to the study of oxidation of Al. The main difference between ReaxFF and the ES+ model according to the published models is the inclusion of dispersi on interactions and the over coordination penalty functions in ReaxFF. Both of these interactions are considered in this work. The coordination functions are set to only affect over coordinated anions in the oxide, which, in turn, destabilizes the NaCl pha se. The best fit for the dispersion interactions is zero in this work, which suggests any dispersion contributions are already incorporated in other interactions. In this work, the Bixby ite phase is penalized with angular terms that also penalize all other phase. In future work, a better solution may arise with a bit more emphasis on the other oxide phases and a better representation of dispersion interactions.
167 A B Figure 51. The crystal structure of Al2O3 (A) has a hexagonal unit cell of spac e group c R 3 Most empirical potentials incorrectly predict the cubic Bixbyite (B) crystal structure with space group 3 Ia to be the ground state structure of Al2O3
168 60 80 100 120 140 160 Bond Angle (Degrees) Bixbyite CorundumAl-O-Al Bond Angle Distribution Figure 5 2 The Al O Al bond angle distribution i n corundum (blue) and Bixbyite (orange).The bond angles around O in the Bixbyite structure are closer to the ideal tetrahedral angle of 109.47o than in corundum.
169 Figure 5 3 The energy vs. volume of several aluminum oxide phases. The ground state is the phase (black). Solid black point are determined with DFT (PBE). Open symbols are values determined using the COMB potential Bixbyite ( b lue) is shifted above the ground state but lies below Al2O3 (green). AlO in the NaCl lattice (red) is shifted up with a coordination penalty function.
170 Figure 54 A closed packed plane of Al atoms can add to the Al2O3 (0001) at three locations. Relative to the oxygen sublattice, atoms can be added with an FCC (F) orientation, HCP (H) orientation or directl y on top (O) of the O atoms on the surface.
171 A. H OAl Al2 1 B. F OAl Al 2 1 Figure 5 5 A second layer of Al atoms in the FCC metal film can be oriented over either of the two remaining sites. The surface in A is formed by aligning the first Al layer in the metal over the O sites,OAl1 and the second Al layer over the H sites, HAl2. The surface in B has the second Al layer aligned over the F sites, FAl2.
172 A. O HAl Al2 1 B. F OAl Al2 1 Figure 5 6 The interfaces with the highest (A) and lowest (B) W sep as predicted by the COMB potential Color corresponds to charge ranging from red=1.5 and blue= 1.5 formal charge units.
173 -20 -15 -10 -5 0 5 10 15 20 -1 0 1 2 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 Planar Distribution FunctionDistance () Averaged Charge Density per 2 Figure 5 7. The charge density (blue line) and separation distance per plane of atoms across the O HAl Al2 1 interface at 0 K. The red line represents the distance between planes of O atoms and black corresponds to Al.
174 Table 5 1. Atomic and e lectrostatic p otential p arameters Al O (eVq 1 ) 4.478988 4.700782 J (eVq 2 ) 3.246069 5.064537 K (eVq 3 ) 0.000000 2.756183 L (eVq 4 ) 0.000000 0.992188 ( 1 ) 0.618852 3.012029 (q) 0.146428 0.030819 P ( 3 ) 0.335000 0.323757 JP1(eVq3r5) 0.612118 0.054039 JP1( eVq4r5) 0.006107 1.136518 D u () 0.011615 1.628749 D l () 0.679563 0.244020 Q u (q) 3.000000 6.000000 Q l (q) 5 .00000 0 2.000000 CN* 6 .00000 0 4. 0 00000 P VdW 0 .00000 0 0 0 00000 Coord 0.500000 0. 250000 E Coord 0.000000 4.887383 Table 52. Interaction d ependent p otential p arameters Al Al O O Al O O Al A (eV) 379.6596 3 523.35 9 1023.282 1023.282 B (eV) 70.79098 204.625 9 132.1803 132.1803 ( 1 ) 2.197639 5.516839 1.399757 1.39 9757 ( 1 ) 1.223328 2.527568 3.452255 3.452255 0.168176 2.00000 0 0.168176 2.00000 0 1.000000 1.00000 0 1.000000 1.00000 0 M 1.000000 1.00000 0 1.000000 1.00000 0 C (rad.) 0.000705 43.56000 1.832422 0.000747 D (rad.) 0.305736 1.000000 1.000000 1.0000 0 H (rad.) 0.047803 0.220000 0.131031 0. 415250 n B 10.00000 10.0000 0 10.00000 10.00000 Rs () 3.000000 3.0 000 00 2.200000 2.200000 Ss () 3.400000 3.4 000 00 2.600000 2.600000
175 Table 53. Coefficients for Legendre p olynomials Al Al Al Al O Al O Al O 1 LPK (eV) 0.015581 0.165976 0.477644 2 LPK (eV) 0.039822 0.279429 0.966470 3 LPK (eV) 0.072865 0.247754 0.050423
176 Table 54. Properties of m etallic aluminum Experiment or DFT fES+ gReaxFF COMB FCC Properties a ( 4.05 a 4.0 5 4.01 4.0 5 E o (eV/atom) 3.36 b 3. 39 3. 36 3. 36 B (GPa) 79 c 83 79 79 C 11 (GPa) 114 c 94 119 114 C 12 (GPa) 62 c 77 57 62 C 44 (GPa) 32 c 34 50 36 Phase Transitions (eV/atom) E o HCP 0.03 d 0.00 0.03 E o BCC 0. 11 d 0.05 0.12 E o Diamond 0.67 d 0.70 1.38 E o Cubic 0.33 d 0.37 0.61 H f Point Defects (eV) V Al 0.68 d 0.85 1.52 I oct 2 .79 d 3.33 100 dumbell 2.50 110 dumbell 7.69 111 dumbell 8.58 Planar Defects (mJm 2 ) (100 ) 980 e 943 g 576 948 (110) 980 e 1000 g 576 1100 (111) 980 e 870 g 576 805 ISF 166 e ,120 144 h 135 USF 168 d 858 Twin 75 e 68 Indicates predicted values not included in the training set. aReference 157 bReference 143 cReference 172 dReference 140 eReference 196 fReference 22 gReference 40 hReference 197, 198
177 Table 5 5. Ionizat ion p otentials and e lectron a ffinities (in eV) for Al and O Al Ab Initio a COMB EA 1. 23 1.08 1 st IP 7. 73 7. 72 2 nd IP 12.6 14.2 3 rd IP 36.6 20.7 O 1 st IP 13.1 13.1 1 st EA 1.40 1.40 2 nd EA 6.08 6.08 aCalculated with CCSD(t)/cc pVTZ Table 56. P r operties of Al2O3 Experiment 1 st Principles fES+ gReaxFF COMB a ( 4. 76 a 4.76 4.81 4. 75 c ( 13.0 b 13.0 13.1 12.9 E coh (eV/ Al 2 O 3 ) 31.8 b 31.8 38 31.9 Bulk Mod. (GPa) 254 c 250 248 256 C 11 (GPa) 497 c 537 546 C 12 (GPa) 164 c 180 175 C 13 (GPa) 111 c 106 88.2 C 1 4 (GPa) 24 c 30 27.8 C 33 (GPa) 498 c 509 507 C 44 (GPa) 147 c 130 195 C 66 (GPa) 168 c 179 201 Al Charge 2.86 0.95 Defects in Al 2 O 3 (eV) V Al 8.44 d 4.46 V O 6.09 d 0.90 O i 8.84 d 3.61 Al i 16.8 d 10.8 (0001) (J m 2 ) 1.5 9/2.12 e 2.67 1.0 1.82 E for Meta Stable Oxides (eV/Al 2 O 3 ) Bixbyite 0.97 a 0.66 0.74 0.38 a 0.19 1.84 aReference 199 bReference 143 cReference 200 dReference 194 eReference 184 (GGA/LDA) fReference 22 gReference 40
178 Table 57. Interfacial separation distance, Do (in ) and Wsep (in J.m2) for the relaxed single Al terminated Al(111)/ Al2O3(0001) interfaces. O FAl Al2 1 H FAl Al2 1 F HAl Al2 1 O HAl Al2 1 F OAl Al2 1 H OAl Al2 1 a W sep Expt 1. 13 b W sep DFT* 1. 36 1. 36 0. 69 0. 69 1.18 1.18 c W sep ReaxFF* 1.74 1.74 1.31 1.31 W sep COMB 1.31 1.31 1.31 1.39 1.29 1. 38 b D o DFT* 0.70 0.70 2.57 2.57 1.62 1.62 D o COMB 2.55 3.01 2.55 2.54 2.56 2.54 *Values do not distinguish between Al2 orientations aReference 86 bReference 184 cReference 40
179 CHAPTER 6 CURRENT CAPABILITIES AND FUTURE POSSIBILI TIES OF CHARGE OPTIMZED MANY BODY POTENTIALS Development of empirical potentials is an active area of research in computational materials science and engineering. This is particularly true since device design has reached the nanoscale where fully atomistic simulations of actual life sized systems are feasible. Often such applications involve the interface between two different functional materials. As the size of a multifunctional device is reduced, the material contained within the surface and interfacial regions const itutes an increased portion of the total device and, commensurately, performance is increasingly influenced by surface and interfacial effects.173 Computational simulations hold promise for providing insight into complex interfacial and surface processes and their impact on device performance. Computational applications s uch as these rely heavily upon classical analytical potentials. The availability of suitable potential energy functions is currently the main limitation to simulations of realistically complex systems Chapter 3 of this work presents a comparison between the empirical REBO potential and first principle methods. The primary conclusion to draw, in regards to methodology, is all levels of theory have significant strengths and weakness that must be balanced and accounted for. In that work, the higher quality m ethods are limited to small systems and low quality statistics, yet are more reliable in their predictions. DFT results suggest that CF3 and CF3+ react differently at deposition energies of 50 eV. At that energy, CF3 adds to the surface and forms bonds w hile CF3+ fragments and removes more H from the surface. The REBO potential does not include electrostatic interactions and so can not resolve the difference in reactivity between these species.
180 An assumption is that at high deposition energies, both speci es will react similarly. This assumption is shown to be invalid at 50 eV. DFT MD also reveals a single step addition of CF02 to the diamond surface when coupled with the formation and evolution of HF. The REBO potential failed to predict this key react ion and several other HF forming reactions. The discrepancy is attributable to the short interaction distance between H and F in REBO. The bond lengths in key transition states are found to be at the extreme limits of interaction distance in REBO. The efficiency of REBO allows for a sampling size that is several orders of magnitude larger than DFT methods. The improved statistics with the REBO model have the potential to significantly enhance the quality of the results. Unfortunately, through compari son to higher level methods, the limitations of the model and their implications in the results are revealed. More importantly, weaknesses in the model can be addressed in the next generation of the methodology, which for the REBO potential, involves the addition of variable charge electrostatics and the increased interaction distances. The extension of COMB to hydrocarbons is currently an active project. As the COMB model is extended to systems with challenging chemistry and an abundance of structures, the transferability of the method will be challenged. The current version of the model is designed to give the best chance of success in expanding to these systems. As shown in Chapter 2, the potential is derived from the same DFT functionals as the current hydrocarbon potentials.16, 24, 25 Brenner attributes the transferability of REBO to this basis in DFT.52 Furthermore, additional functions that, so
181 far, have been required to fit new systems, are chosen to correct interactions that are lost in the approximations of the initial derivation. This project initially began with three potentials, Yasukawas potential 17, Streitz and Mintmires ES+ model ,22 and a modified embedded atom method (MEAM) .201, 202 The idea at the time was to compare these three potentials for application to multi phase systems. The Yasukawa potential is essentially a Tersoff potential coupled with a variable charge scheme. The short range interactions vary with charge. However, in that potential, the electrostatic functions are based on a point charge model with interactions cutoff just past the second neighbor shell, which causes several instabilities in application. The ES+ model has electrostatic interactions are determined as functionals of charge density with truly long interaction distances. Unfortunately, the electrostatics are couples to and Finnis Sinclair short range potential which cannot describe covalent systems. MEAM is an improved second moment model that includes angular terms to model covalent systems but lacks electrostatic interactions. 201 Of the three approaches, the Yasukawa potential, while flawed, performed the most reliably. With the hindsight gained from the derivation in Chapter 2, it is clear that the Yasukawa formalism succeeds because it adds back some of the interactions that are approximated away in ES+ and MEAM. ES+ is based on the rigid atom approximation and, as a consequence, is unreliable in complex systems. Yasukawa addresses the short range effects that are neglected with the rigid ion approximation. In the COMB formalism, the long range electrostatic interactions that are missing in Yasukawas potential 17 are added back along with several long range i nteractions that are lost with
182 a rigid atom approximation. After much trial in error, it seems this current course of development is the most promising. The adaptability of the COMB potential is demonstrated by fitting the potential for two metal metal oxide systems. First, demonstrated in Chapter 4 is the development of a potential for the Cu/CuO system, which involves multiple oxidation states of Cu. The initial oxide phase formed under oxygen lean conditions is the metastable cubic phase of Cu2O. As oxidation proceeds and at higher temperatures, monoclinic CuO forms. The potential is parameterized to model the initial stages of oxidation through the growth of Cu2O islands. The results of the fit suggest the potential is well suited for this task. The bulk and mechanical properties are accurately reproduced compared with experimental and first principles based data. Surface energies and defect formation enthalpies are also consistent with first principles calculations. The potential reproduces the relative phase order and formation enthalpies for the low energy oxide phases. The potential captures the dissociation energy of molecular oxygen and its anions, which, in turn, leads to a prediction of adsorption energy consistent with DFT. The potenti al also predicts an adsorption induced missing row reconstruction as the stable oxidized surface at standard temperature, consistent with experiment. The ability of this potential formalism to be adapted to both metallic phases and the oxides both as soli d phases and single molecules and with reconstructed surfaces demonstrates a remarkable flexibility in the functional form. In Chapter 5, the potential is adapted to model Al/ Al2O3 interfaces. The low energy phase of the oxide, Al2O3, forms in the hexagonal corundum crystal structure, which is stabilized by subtle polarization effects in the electron density distribution. Most
183 classical potentials fail to predict the correc t ground state crystal structure and instead predict a cubic Bixbyite structure. The COMB model is able to correctly predict corundum as the ground state structure. This is a distinction that is not easily made even with first principles methods. Similar to Cu/CuO system, the potential also reproduces the bulk and mechanical properties of both metallic and oxide phases as well as the surface energies and defect formation enthalpies in both phases. It predicts atomically sharp interfaces between Al(111) surface and crystalline Al2O3(0001), which is consistent with the results of DFT calculations. However, the potential does not reproduce interfacial reconstructions as predicted by DFT. Interfacial works of separation are comparable to DFT calculated val ues. However, trends in relative energies of different variations of the interface are not captured, indicating the potential is best used for bulk phase calculations and perhaps studies of oxidation in this material. In other efforts the potential has been extended to model all the low energy polymorphic phases of SiO218 including amorphous phases and Hf HfO2184. These systems also pose unique challenges that have not been addressed with less flexible models. For example, most potentials for HfO2 fail to predict the correct phase order.184 The flexibility and transferability of the COMB formalism as demonstrated here and in these other works is quite encouraging. 6.1 The COMB Code The COMB project was initiated to address the growing need for n ew efficient computational methods. The implied tasks in that objective are to first develop a potential formalism with sufficient flexibility to model many materials of interest and secondly to concurrently develop the tools and procedures to efficiently fit the potential
184 to new materials. The two tasks are somewhat interrelated where now the COMB project consists of both a potential form and a suite of FORTRAN90 codes. In addition to the flexibility of the current model, the code itself is designed with the machinery in place to efficiently add additional interaction as the needs arise. For example, a possible improvement for Al2O3 would be to add the quadrapole moments. The code to add higher order polarization modes along the lines of the compressible ion models26 has already been added with the point dipole model. There are some technical issues to address, such as how polarizability varies with charge. However, the administrative coding of how the interactions are processed and information passed is complete. For entirely new interactions such as four bodied energy contributions that arise in hydrocarbon models, the architecture to pass information, even in parallel is in place. The loop structure would have to be created for the four bodied interactions. However, the code is modularized and fully declared with no implicit variables to assist modification The current code exists as two versions: one parallel and one serial. The parallel version breaks the calculation up into equal spatial regions. For example, a 2x2 spatial decomposition breaks the system up in two dimensions into a 2x2 grid. Each section of the grid is assigned a processor that only processes that part of the system. A 2x2 system should, therefore, process nearly four times faster than the same system on one processor. The parallel code is designed for running large systems with the focus on computational efficiency. The code scales linearly with system size depending on the number of processor used.
185 A serial code processes the calculation on one processor, which allows for smaller systems that are used mainly for potential development. Fitting a potential requires many sequential calculations. A well designed serial code efficiently processes each calculation as quickly as possible, where as, in a parallel code; emphasis is placed on processing many calculations simultaneously. The two code formats together are able to handle both development and production. 6.2 The Fitting Problem The models developed here are fit to a mix of measured and calculated values. The preference is to fit to empiric al data augmented with first principles calculations when empirical values are not obtainable. This is a philosophical choice. Several potentials are fit entirely to calculated data sets.16, 202 However, for the systems modeled here, the first principles values are not always accurate or correct.144 There are arguments for either approach.140 A model fit to experimental data is fit to a more accurate data s et than what is achievable with first principles calculations. A model fit with precision to the more accurate data set will reproduce those values more accurately. However, experimental data sets are small and transferability beyond the data set is never ensured, even with a precise fit. In contrast, the only limit on the size of a training date set based on first principles calculations is computer time. Given enough resources, a training data set that maps the DFT potential energy surface could, in princ iple, be generated. Any model fit to such an extensive data set could be as transferable as DFT, yet only as accurate as DFT. With the current state of the art, the best approach is to fit to a mix of experimental and calculated data points. Brenner concl uded in his review of bond order potentials that while it is quite possible that no purely empirical potential will ever be truly transferable, no potential derived entirely
186 from first principles will be truly accurate.52 Brenners point is appreciated in this work considering systems such as CuO, where being as accurate as DFT is not very accurate. Perhaps, as important as developing the potential is designing the fitting process. A large part of the success of the ReaxFF models is that it has been extended to s o many systems. Within the COMB suite are series of codes designed to do multi dimensional minimizations. The algorithms are standards from the numerical recipes.203 The fitting codes work with the actual serial version of the COMB code. They call the same energy and force subroutines that are used in the production codes and use the same variable names and data types so that parameters can be quickly tested and put to use. Each phase of a multiphase system can be fit separately. For example all metal metal interactions can be fit to bulk metallic phases. Self energies parameters are fit to atomic values. The different interaction can also be optimized separately; charge dependent and charge independent parameters can be optimized to a degree separately. Versions of the fitting codes have been optimized to specifically fit these separate phases. Fitting a potential is a challenging task. The COMB suite has evolved to where the basic tools need to run simulations, analyze results and fit new systems are in place. 6.3 Future Applications In addition to hydrocarbon systems, the COMB model is currently being extended to Zn and ZnO, U and UO2, Ti and TiO4. This gives an idea of the diverse range of materials and applications to which t he model can be applied. The next phase of development will be to extensively run applications and comparative studies. This is an
187 area with fairly sparse data, as the model has just recently been finalized. The best test is to get the code into producti on and get feedback on its strengths and failings. As it stands, the COMB model and suite has the potential to greatly expand computation to systems that cannot be studied by any other means, which in turn, could tremendously enhance the choice and select ion of materials for the next generation of multi functional devices.
188 LIST OF REFERENCES 1 P. Mutti, G. Ghislotti, S. Bertoni, et al., Appl. Phys. Lett. 66 851 (1995). 2 T. Shimizu Iwayama, S. Nakao, and K. Saitoh, Appl. Phys. Let t. 65 1814 (1994). 3 Q. Zhang, S. C. Bayliss, and D. A. Hutt, Appl. Phys. Lett. 66, 1977 (1995). 4 L. Hanley, Y. Ch oi, E. R. Fuoco, et al. Nucl, Inst. Meth Phys Res B 203, 116 (2003). 5 K. Kimura and S. Iwasaki, J. Appl. Phys. 83, 1345 (1998). 6 P. F. Trwoga, A. J. Kenyon, and C. W. Pitt, J. Appl. Phys. 83, 3789 (1998). 7 K. S. Seol, A. Ieki, Y. Ohki, et al., J. Appl. Phys. 79, 412 (1996). 8 M. L. Brongersma, A. Polman, K. S. Min, et al. Appl. Phys. Lett. 72 2577 (1998). 9 K. S. Zhuravlev, A. M. Gilinsky, and A. Y. Kobitsky, Appl. Phys. Lett. 73, 2962 (1998). 10 R. A. Buckingham, Proc. R. Soc. Lond. A 168 264 (1938). 11 J. Tersoff, Phys. Rev. B 37, 6991 (1988). 12 J. Tersoff, Phys. Rev. B 38, 9902 (1988). 13 B. J. Garrison, E. J. Dawnkas ki, D. Srivastava, et al. Science 255, 835 (1992). 14 A. K. Rappe and W. A. Goddard, J. Phys. Chem. 95, 3358 (1991). 15 T. Iwasaki and H. Miura, J. Mater. Res. 16, 1789 (2001). 16 A. C. T. van Duin, S. Dasgupta, F. Lorant, et al. J. Phys. Chem. A 105, 9396 (2001). 17 A. Yasukawa, J SME Intl. J. A 39, 313 (1996). 18 M. B. J. Wijesundara, Y. Ji, B. Ni, et al. J. Appl. Phys. 88, 5004 (2000). 19 Q. Zhang, T. Cagin, A. van Duin, et al. Phys. Rev. B 69 (2004). 20 X. W. Zhou and F. P. Doty, Phys. Rev. B 78 (2008). 21 M. W. Finnis and J. E. Sinclair, Philo. Mag. A 50, 45 (1984).
189 22 F. H. Streitz and J. W. Mintmire, Phys. Rev. B 50, 11996 (1994). 23 G. C. Abell, Phys. Rev. B 31, 6184 (1985). 24 D. W. Brenner, Phys. Rev. B 42, 9458 (1990). 25 D. W. Bren ner, O. A. Shenderova, J. A. Harrison, et al. J. Phy s. : Conden. Matter 14, 783 (2002). 26 T. R. Shan, B. D. Devine, T. W. Kemper, et al. Phys. Rev. B 81 (2010). 27 W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 A1138 (1965). 28 P. Hohenberg and W. Kohn, Phys. Rev. B 136 B864 (1964). 29 A. J. Cohen, P. Mori Sanchez, and W. T. Yang, Science 321, 792 (2008). 30 R. G. Parr and W. Yang, Density Functional Theory of Atoms and M olecules (Oxford University Press, New York, 1989). 31 R. H. Milne and A. Howie Philo. Mag. A 49, 665 (1984). 32 R. Guan, H. Hashimoto, and T. Yoshida, Acta Crystal. B 40, 109 (1984). 33 L. O. Brockway and I. M. Adler, J. Electrochem. Soc. 119, 899 (1972). 34 A. Ronnquist and H. Fischmeister, J Inst. of Metals 89, 65 (1960). 35 R. Drautz, D. A. Murdick, D. Nguyen Manh, et al. Phys. Rev. B 72 (2005). 36 D. G. Pettifor, M. W. Finnis, D. Nguyen Manh, et al ., Mate r. Sci. Eng. A 365, 2 (2004). 37 G. W Zhou and J. C. Yang, Mater. at High Temperatures 20, 247 (2003). 38 M. Schvar tzman, A. Math ur, Y. Kang, et al ., J. Vac. Sci. Tech. B 26, 2394 (2008). 39 B. Devine, I. Jang, T. Kemper, et al ., J. Phys. Chem. C 114, 12535 (2010). 40 Q. Zhang, T. Cagin, A. van Duin, et al ., Phys. Rev. B 69 (2004). 41 M. Wilson, M. Exner, Y. M. Huan g, et al. Phys. Rev. B 54, 15683 (1996). 42 D. G. Pettifor, Phys. Rev. Lett. 63, 2480 (1989). 43 Oleinik, II and D. G. Pettifor, Phys. Rev. B 59 8500 (1999).
190 44 D. Raymand, A. C. T. van Duin, M. Baudin, et al ., Surf. Sci. 602, 1020 (2008). 45 D. A. M urdick, X. W. Zhou, H. N. G. Wadley, et al ., Phys. Rev. B 73 (2006). 46 B. A. Gillespie, X. W. Zhou, D. A. Murdick, et al ., Phys. Rev. B 75 (2007). 47 R. Drautz, X. W. Zhou, D. A. Murdick, et al ., Prog. Mater. Sci. 52, 196 (2007). 48 S. Znam, D. NguyenManh, D. G. Pettifor, et al ., Philo. Mag. 83, 415 (2003). 49 S. Ryu, C. R. Weinberger, M. I. Baskes, et al. Modelling and Simulation in Materials Science and Engineering 17 (2009). 50 J. Uddin, M. I. Baskes, S. G. Srinivasan, et al ., Phys. Rev. B 81 (20 10). 51 T. Hammerschmidt, R. Drautz, and D. G. Pettifor, Intl. J. Mater. Res. 100, 1479 (2009). 52 D. W. Brenner, Phy s. Status Solidi B 217, 23 (2000). 53 P. A. M. Dirac, Proc. R. Soc. Lond. A 117, 610 (1928). 54 A. Szabo and N. Ostlund, Modern Quantum C hemistry (Macmillan, New York, 1982). 55 P. A. M. Dirac, Proc. R. Soc. Lond. A 112, 661 (1926). 56 W. Heisenberg, Zeitschrift Fur Physik 40, 501 (1926). 57 M. Born and R. Oppenheimer, Annalen Der Physik 84, 0457 (1927). 58 R. C. Brown, C. J. Cramer, and J. T. Roberts, J. Phys. Chem. B 101, 9574 (1997). 59 M. Born and K.Huang, Dynamical Theory of Crystal Lattices (Oxford University Press, Oxford, 1954). 60 Ditchfie.R, W. J. Hehre, and J. A. Pople, J. Chem. Phys. 54, 724 (1971). 61 P. A. M. Dirac, Pr oc. R. Soc. Lond. A 111, 405 (1926). 62 W. Pauli, Zeitschrift Fur Physik 31, 765 (1925). 63 J. C. Slater, Phys. Rev. 35, 0210 (1930).
191 64 P. A. M. Dirac, Proc. Camb. Philo. Soc 26, 376 (1930). 65 D. R. Hartree, Proc. Camb. Philo. Soc 24, 89 (1928). 66 V. Fock, Z. Phys. 61, 126 (1930). 67 R. S. Mulliken, J. Chem. Phys. 23, 2338 (1955). 68 R. S. Mulliken, J. Chem. Phys. 23, 2343 (1955). 69 R. S. Mulliken, J. Chem. Phys. 23, 1833 (1955). 70 R. S. Mulliken, J. Chem. Phys. 23, 1841 (1955). 71 P. O. Lowd in, Adv. Quantum Chem. 5 185 (1970). 72 A. Freedman and C. D. Stinespring, Appl. Phys. Lett. 57 1194 (1990). 73 C. A. Coulson, Proc. R. Soc. Lond. A 169, 413 (1939). 74 P. E. Blochl, Phys. Rev. B 50, 17953 (1994). 75 E. Fermi, Atti Accad. Naz. Lincei Cl. Sci. Fis., Mat. Nat., Rend. 6 602 (1927). 76 A. K. George, A. S. Harry, B. J. Berne, et al., J. Comp. Chem. 23, 1515 (2002). 77 R. P. Feynman, Phys. Rev. 56, 340 (1939). 78 W. Matthew, C. Foulkes, and R. Haydock, Phys. Rev. B 39, 12520 (1989). 79 J. Harris, Phys. Rev. B 31, 1770 (1985). 80 B. G. Dick and A. W. Overhauser, Phys. Rev. 112, 90 (1958). 81 D. Wolf, P. Keblinski, S. R. Phillpot, et al., J. Chem. Phys. 110, 8254 (1999). 82 S. W. Rick, S. J. Stuart, and B. J. Berne, J. Chem. Phys. 101, 6141 (1994). 83 W. C. Swope, H. C. Andersen, P. H. Berens, et al ., J. Chem. Phys. 76, 637 (1982). 84 F. Cyrot Lackmann, J. Phys. Chem. Solids 29, 1235 (1968). 85 D. W. Brenner, Phys. Rev. Lett. 63, 1022 (1989). 86 I. Batyrev, A. Alavi, and M. W. Finn is, Faraday Discussions, 33 (1999).
192 87 L. A. Curtiss, K. Raghavachari, P. C. Redfern, et al ., Chem. Phys. Lett. 314, 101 (1999). 88 G. D. Purvis and R. J. Bartlett, J. Chem. Phys. 76, 1910 (1982). 89 M. Sternberg, P. Zapol, and L. A. Curtiss, Molec. Phy s. 103, 1017 (2005). 90 B. T. Thole, Chem. Phys. 59, 341 (1981). 91 J. Applequi st J. R. Carl, and K. K. Fung, J. Am. Chem. Soc. 94, 2952 (1972). 92 D. Va nbelle and S. J. Wodak, Comp. Phys. Comm. 91, 253 (1995). 93 P. Geerlings, F. De Proft, and W. Langenaeker, Chem. Rev. 103, 1793 (2003). 94 H. Toufar, K. Nulens, G. O. A. Janssens, e t al ., J. Phys. Chem. 100 15383 (1996). 95 T. H. Dunning, J. Chem. Phys. 90, 1007 (1989). 96 J. M. Soler, E. Artacho, J. D. Gale, et al. J. Phys: Conden. Matter 14 2745 (2002). 97 E. Artacho, E. Anglada, O. Dieguez, et al ., J. Phys.: Conden. Matter 20, 064208 (2008). 98 I. Jang and S. B. Sinnott, J. Phys. Chem. B 108, 18993 (2004). 99 L. Hanley and S. B. Sinnott, Surf. Sci. 500, 500 (2002). 100 S. A. Adelman and J. D. Doll, J. Chem. Phys. 64, 2375 (1976). 101 R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985). 102 W. G. Hoover, Phys. Rev. A 31, 1695 (1985). 103 S. Nose, J. Chem. Phys. 81, 511 (1984). 104 S. Nose, Molec. Phys. 52, 255 (1984). 105 N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 (1991). 106 J. A. Pople, M. HeadG ordon, and K. Raghavachari, J. Chem. Phys. 87, 5968 (1987).
193 107 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996). 108 A. G. Baboul, L. A. Curtiss, P. C. Redfern, et al ., J. Chem. Phys. 110, 7650 (1999). 109 G. W. T. M. J. Frisch, H. B. Schlegel, G. E. Scuseria, et al (Gaussian, Inc., Wallingford CT., 2004). 110 A. D. Becke, J. Chem. Phys. 98, 5648 (1993). 111 Y. Zhao and D. G. Truhlar, J. Phys Chem. A 112, 1095 (2008). 112 A. D. McL ean and G. S. Chandler, J. Chem. Phys. 72, 5639 (1980). 113 A. D. Boese and J. M. L. Martin, J. Chem. Phys. 121, 3405 (2004). 114 E. I. Izgorodina, D. R. B. Brittain, J. L. Hodgson, et al. J. Phys. Chem. A 111, 10754 (2007). 115 H. L. Schmider and A. D. Becke, J. Chem. Phys. 108, 9624 (1998). 116 K. E. Riley, B. T. Holt, and K. M. Merz, J. Chem Theory Comp. 3 407 (2007). 117 D. G. Goodwin, Appl. Phys. Lett. 59, 277 (1991). 118 S. J. Harris, Appl. Phys. Lett 56, 2298 (1990). 119 M. Page and D. W. Brenner, J. Am. Chem. Soc. 113, 3270 (1991). 120 D. G. Goodwin, J. Appl. Phys. 74, 6895 (1993). 121 S. F. Boys and F. Bernardi, Molec. Phys. 19 553 (1970). 122 L. Pastewka, P. Pou, R. Perez, et al. Phys. Rev. B 78 (2008). 123 M. Mrovec, M. Moseler, C. Elsasser, et al ., Prog. Mater. Sci. 52, 230 (2007). 124 F. W. Young Jr, J. V. Cathcart, and A T. Gwathmey, Acta Metal. 4 145 (1956). 125 K. Heinemann, D. B. Rao, and D. L. Douglass, Oxidation of Metals 9 379 (1975). 126 G. W. Zhou and J. C. Yang, Surf. Sci. 531, 359 (2003). 127 N. Cabrera and N. F. Mott, Reports on Progress in Physics 12, 163 (1948). 128 R. H. Milne and A. Howie, Philo. Mag. A 49, 665 (1984).
194 129 K. Tanaka, T. Inomata, and M. Kogoma, Thin Solid Films 386, 217 (2001). 130 L. Sun and J. C. Yang, J. Mater. Res. 20, 1910 (2005). 131 R. H. Milne, Surf. Sci. 232, 17 (1990). 132 G. Zhou, J. Appl. Phys. 105, 104302 (2009). 133 K. Chenoweth, A. C. T. van Duin, P. Persson, et al ., J. Phys. Chem. C 112, 14645 (2008). 134 R. W. G. Wyckoff, Crystal S tructures (Interscience Publishers, New York, 1963). 135 N. E. Brese, M. Okeeffe, B. L. Ramakrishna, et al ., J. Solid State Chem. 89, 184 (1990). 136 K. Takahashi, A. Itoh, T. Nakamura, et al ., Thin S olid Films 374, 303 (2000). 137 E. Ruiz, S. Alvarez, P. Alemany, et al. Phys. Rev. B 56, 7189 (1997). 138 G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993). 139 J. H. Rose, J. R. Smith, F. Guinea, et al ., P hys. Rev. B 29, 2963 (1984). 140 Y. Mishin D. Farkas, M. J. Mehl, et al ., Phys. Rev. B 59 3393 (1999). 141 A. Soon, X. Y. Cui, B. Delley, et al ., Phys. Rev. B 79 (2009). 142 H. Raebiger, S. Lany, and A. Zunger, Phys. Rev. B 76 (2007). 143 D. R. Lide, C RC Handbook of Chemistry and P hysics (CRC Press, Boca Raton, FL, 2005). 144 D. O. Scanlon, B. J. Morgan, and G. W. Watson, J. Chem. Phys. 131 (2009). 145 M. Sprik, J. Phys. Chem. 95, 2283 (1991). 146 M. Nolan, Thin Solid Films 516, 8130 (2008). 147 A. Soon, M. Todorova, B. Delley, et al ., Phy s. Rev. B 75 (2007). 148 S. M. Woodley, P. D. Battle, C. R. A. Catlow, et al ., J. Phys. Chem. B 105, 6824 (2001).
195 149 K. J. Blobaum, D. Van Heerden, A. J. Wagner, et al. J. Mater. Res. 18, 1535 (2003). 150 F. Besenbacher and J. K. Norskov, Prog. Surf. Sci. 44, 5 (1993). 151 E. A. Colbourn and J. E. Inglesfield, Phys. Rev. Lett. 66, 2006 (1991). 152 T. Fujita, Y. Okawa, Y. Matsumoto, et al., Phys. Rev. B 54, 2167 (1996). 153 H. Iddir, D. D. Fong, P. Zapol, et al ., Phys. Rev. B 76 (2007). 154 K. W. Ja cobsen and J. K. Norskov, Phys. Rev. Lett. 65, 1788 (1990). 155 F. Jensen, F. Besenbacher, E. Laegsgaard, et al ., P hys. Rev. B 42 9206 (1990). 156 T. Kangas and K. Laasonen, Surf. Sci. 602, 3239 (2008). 157 M. Kittel, M. Polcik, R. Terborg, et al Sur f. Sci. 470, 311 (2001). 158 K. Lahtonen, M. Hirsimki, M. Lampimki, et al. J. Chem. Phys. 129, 124703 (2008). 159 T. Lederer, D. Arvanitis, G. Comelli, et al ., P hys. Rev. B 48, 15390 (1993). 160 G. H. Liu and R. G. Parr, J. Am. Chem. Soc. 117, 3179 (1995). 161 I. Merrick, J. E. Inglesfield, and H. Ishida, Surf. Sci. 551, 158 (2004). 162 M. Ahonen, M. Hirsimaki, A. Puisto, et al ., Chem. Phys. Lett. 456, 211 (2008). 163 I. K. Robinson, E. Vlieg, and S. Ferrer, Phys. Rev. B 42 6954 (1990). 164 M. Sotto, Surf. Sci. 260, 235 (1992). 165 S. Stolbov, A. Kara, and T. S. Rahman, Phys. Rev. B 66 (2002). 166 S. Stolbov and T. S. Rahman, Phys. Rev. Lett. 89 (2002). 167 M. Wuttig, R. Franchy, and H. Ibach, Surf. Sci. 224, L979 (1989). 168 M. Wuttig, R. F ranchy, and H. Ibach, Surf. Sci. 213, 103 (1989). 169 H. C. Zeng, R. A. McFarlane, and K. A. R. Mitchell, Surf. Sci. 208, L7 (1989). 170 X. Duan, O. Warschkow, A. Soon, et al ., Phys. Rev. B 81 (2010).
196 171 J. K. Barton, A. A. Vertegel, E. W. Bohannan, et al. Chem Mater. 13 952 (2001). 172 G. Simons and H. Wang, Single Crystal Elastic Constants and Calculated Aggregate P roperties (MIT Press, Cambridge, MA, 1977). 173 R. W. Balluffi and A. P. Sutton, Intergranular and Interphase Boundaries in Materials Part 1, 207 1 (1996). 174 Y. Mishin, M. J. Mehl, D. A. Papaconstantopoulos, et al. Phys. Rev. B 63, 224106 (2001). 175 M. M. Beg and S. M. Shapiro, Phys. Rev. B 13, 1728 (1976). 176 J. Ha llberg and R. C. Hanson, Phys. Status Solidi 42, 305 (1970). 177 V. E. Henrich and P. A. Cox, The Surface Science of Metal O xides (Cambridge University Press, New York, 1994). 178 N. P. Padture, M. Gell, and E. H. Jordan, Science 296, 280 (2002). 179 B. C. Gates, Mate r. Chem. 245, 301 (1995). 180 C. T. Campbell, Surf. Sci. Rep. 27, 1 (1997). 181 J. A. Kelber, Surf. Sci. Rep. 62, 271 (2007). 182 B. Hinnemann and E. A. Carter, J. Phys. Chem. C 111, 7105 (2007). 183 T. J. Campbell, G. Aral, S. Ogata, et al ., Phys. Rev. B 71 (2005). 184 D. J. Siegel, L. G. Hector, and J. B. Adams, Phys. Rev. B 65 (2002). 185 W. Zhang and J. R. Smith, Phys. Rev. Lett. 85, 3225 (2000). 186 J. R. Smith and W. Zhang, Acta Mater. 48, 4395 (2000). 187 T. Campbell, R. K. Kalia, A. Nakano, et al ., Phys. Rev. Lett. 82, 4866 (1999). 188 D. L. Medlin, K. F. McCarty, R. Q. Hwang, et al. Thin Solid Films 299, 110 (1997). 189 G. Paglia, A. L. Rohl, C. E. Buckley et al. J. Mater. Chem. 11, 3310 (2001). 190 G. Paglia, A. L. Rohl, C. E. Buckley, et al. Phys. Rev. B 71 (2005).
197 191 P. Lamparter and R. Kniep, Physica B 234, 405 (1997). 192 M. Vermeersch, R. Sporken, P. Lambin, et al ., Surf. Sci. 235, 5 (1990). 193 U. Schr oder, Solid State Comm. 4 347 (1966). 194 N. D. M. Hine, K. Frensch, W. M. C. Foulkes, et al. Phys. Rev. B 79 (2009). 195 R. Di Felice and J. E. Northrup, Phys. Rev. B 60, R16287 (1999). 196 L. E. Murr, Interfacial Phenomena in Metals and A lloys (AddisonWesley, Reading, MA, 1975). 197 K. H. Westmaco and R. L. Peck, Philo. Mag. 23, 611 (1971). 198 R. H. Rautioaho, Phy s. Status Solidi B 112, 83 (1982). 199 W. H. Gitzen, Alumina as a C eramic M aterial (The American Ceramic Society, Columbus, OH, 1970). 200 J. B. Wachtman, W. E. Tefft, D. G. Lam, et al. J. Res Natl. Bur. Stnds. A. 64, 213 (1960). 201 M. I. Baskes and R. A. Johnson, Modeling and Simulation in Materials Science and Engineering 2 147 (1994). 202 M. I. Baskes, S. G. Srinivasan, S. M. Valone, et al. Phys. Rev. B 75 (2007). 203 W. H. Press, Numerical R ecipes in FORTRAN 77: The A rt of Scientific C omputin g (Cambridge University Press, Cambridge, 1992).
198 BIOGRAPHICAL SKETCH Bryce David Devine received a b achelors degree in chemistry from Boston University in 1992. Following his undergraduate education, he was commissioned as an officer in the US Army and served 5 years as a paratrooper in Ft. Richardson, Alaska. After his military service, Bryce worked as an industrial chemist for several years during which time he was awarded three patents for his work on metal oxide coatings. In the fall 2004, Bryce was afforded the opportunity to pursue his graduate studies at the University of Florida under the advi sement of Prof. Susan Sinnott. He earned his PhD in materials science and engineering in December 2010.