TRANSPORT AND PHASE TRANSITIONS IN CERAMICS AND AMORPHOUS METALS BY ATOMIC LEVEL SIMULATION By BRIAN J. DEMASKE 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 2017
2017 Brian J. Demaske
To m y friends and family
4 ACKNOWLEDGMENTS I want to thank my advisor Prof. Simon Phillpot for his guidance over the past four and a hal f years. He has always been receptive and understanding to any issues that have arisen during my graduate education, both in my research and personal life. He exemplified the role of a mentor by always putting his students and their well being before anything else. His words of encouragement often came at just the right time to give me that much needed boost of confidence. He has promoted me to think independently and gradually pushed me towards self sufficiency in my work. Both skills that are invaluable to a successful career in scientific research. I would also like to acknowledge the other members of my doctoral committee: Prof. Sinnot t Prof. Spearot, Prof. Yang, Prof. Andrew and Prof. Hamlin. I want to thank them for the time and energy they have invested into this process. I would especially like to thank Prof. Spearot for his thoughtful discussions and guidance during early Friday mo rning meetings on shock waves. I would also like to extend my thanks to Prof. Tonks for taking time out of his schedule to serve as a substitute during my oral defense and for helpful suggestions on how to improve this dissertation. Although I am not a mem ber of his group, Prof. Hennig kindly offered me an account on his allocation on the U nive his extra bit of computing power made a world of difference in my work and I thank him for his generosity. I would also like to thank Dr. Alex Chernatynskiy who helped me tackle my first project upon joining the group. Last, but not least, I would like to thank the other students in our group, both past and present, who helped to make graduate school enjoyable. Making it this far would not have been possible without the love and support of my family. My parents gave me the foundation of a strong family and taught me the
5 importance of maintaining a balance between life and work. Both my brother and sister have always been there for me from our e arly years growing up in Ohio all the way up to today. I could not imagine my life without them. Finally, I would like to thank Danielle, who has become my closest friend over these last few years. Her constant love and support have helped me through some rough patches. I could not have done this without her by my side
6 TABLE OF CONTENTS page ACKNOWLEDGMENTS ................................ ................................ ................................ ...... 4 LIST OF TABLES ................................ ................................ ................................ ................ 9 LIST OF FIGURES ................................ ................................ ................................ ............ 11 LIST OF ABBR EVIATIONS ................................ ................................ ............................... 15 ABSTRACT ................................ ................................ ................................ ........................ 16 CHAPTER 1 INTRODUCTION ................................ ................................ ................................ ........ 18 1.1 General Overview ................................ ................................ .............................. 18 1.2 Motivation of Problems ................................ ................................ ...................... 19 1.2.1 Vanadium Carbide Diffusion Barrier ................................ ......................... 19 1.2.2 High strain rate Deformation of Metallic Glasses ................................ .... 21 1.2.3 Lanthanum Borogermanate Glass ceramics ................................ ........... 24 1.3 Outline ................................ ................................ ................................ ................ 25 2 OVERVIEW OF SIMULATION METHODS ................................ ................................ 27 2.1 Density Functional Theory ................................ ................................ ................. 27 2.2 Modern Theory of Polarization ................................ ................................ .......... 29 2.3 Molecular Dynamics ................................ ................................ ........................... 30 2.4 Embedded Atom Method ................................ ................................ ................... 32 3 IN TRINSIC DEFECTS AND SELF DIFFUSION IN VANADIUM CARBIDE ............. 34 3.1 Background ................................ ................................ ................................ ........ 34 3.2 Simulation Methods ................................ ................................ ........................... 36 3.3 Si mulation Results ................................ ................................ ............................. 39 3.3.1 Ordered Structures of V 2 C ................................ ................................ ........ 39 3.3.2 Electronic and Elastic Properties of V 2 C ................................ .................. 44 3.3.3 Vacancy Formation Energies ................................ ................................ ... 48 3.3.4 Energy Barriers for Self diffusion ................................ ............................. 50 184.108.40.206 Vanadium self diffusion ................................ .............................. 50 220.127.116.11 Carbon self diffusion ................................ ................................ ... 56 3.4 Discussion ................................ ................................ ................................ .......... 62 3.5 Summary ................................ ................................ ................................ ............ 68 4 ENERGETICS AND KINETICS OF METAL IMPURITIES IN VANADIUM CARBIDE ................................ ................................ ................................ .................... 69
7 4.1 Background ................................ ................................ ................................ ........ 69 4.2 Simulation Methods ................................ ................................ ........................... 70 4.3 Simulation Results ................................ ................................ ............................. 74 4.3.1 Impurity Site Preference ................................ ................................ ........... 74 4.3.2 Impurity vacancy Binding ................................ ................................ ......... 77 4.3.3 Mi gration Barriers for Impurity Diffusion ................................ ................... 81 4.4 Summary ................................ ................................ ................................ ............ 83 5 SHOCK COMPRESSION OF COPPER ZIRCONIUM METALLIC GLASSES ......... 84 5.1 Background ................................ ................................ ................................ ........ 84 5.2 Simulation Methods ................................ ................................ ........................... 86 5.2.1 Interatomic Potential Verification ................................ .............................. 86 5.2.2 Amorphous Structure Creation ................................ ................................ 89 5.2.3 Piston driven Shock Method ................................ ................................ .... 91 5.3 Simulation Results ................................ ................................ ............................. 93 5.3.1 Shock Wave Structure at Different Shock Intensities .............................. 93 5.3.2 Hugoniot Elastic Limit ................................ ................................ ............... 95 5.3.3 U p U s Shock Hugoniot ................................ ................................ .............. 96 5.3.4 Composition Dependence of the Shock Response ............................... 100 5.4 Discussion ................................ ................................ ................................ ........ 105 5.4.1 Pressure dependent Yield Criterion ................................ ....................... 105 5.4.2 Bilinear U p U s Hugoniot ................................ ................................ .......... 107 5.5 Summary ................................ ................................ ................................ .......... 107 6 MICROSTRUCTURAL RESPONSE OF METALLIC GLASSES TO SHOCK LOADING ................................ ................................ ................................ .................. 109 6.1 Background ................................ ................................ ................................ ...... 109 6.2 Simulation Methods ................................ ................................ ......................... 111 6.3 Simulation Results ................................ ................................ ........................... 114 6.3.1 Shock induced Deformation ................................ ................................ ... 114 6.3.2 Composition Dependence and Melting ................................ .................. 119 6.3.3 Local Structural Changes During Shock Compression ......................... 124 6.4 Summary ................................ ................................ ................................ .......... 129 7 FERROELECTRICTY IN LANTHANUM BOROGERMANATE ............................... 131 7.1 Background ................................ ................................ ................................ ...... 131 7.2 Simulation Methods ................................ ................................ ......................... 132 7.3 Simulation Results ................................ ................................ ........................... 133 7.3.1 Low temperature Phase of LaBGeO 5 ................................ .................... 133 7.3.2 High temperature Phase of LaBGeO 5 ................................ .................... 139 7.3.3 Lattice Dynamics of the Ferroelectric Phase Transition ........................ 142 7.3.4 Spontaneous Polarization ................................ ................................ ....... 146 7.4 Discussion ................................ ................................ ................................ ........ 147 7.5 Summary ................................ ................................ ................................ .......... 149
8 8 CONCLUSIONS ................................ ................................ ................................ ........ 151 LIST OF REFERENCES ................................ ................................ ................................ 158 BIOGRAPHICAL SKETCH ................................ ................................ .............................. 173
9 LIST OF TABLES Table page 3 1 2 C and the low energy structure belonging to space group Pnnm ................................ ... 41 3 2 Zero temperature elastic constants and elastic moduli for the three ordered V 2 C structures. ................................ ................................ ................................ ........ 45 3 3 Vacancy formation ener gies and vibrational formation entropies for the three ordered V 2 C structures calculated using DFT. ................................ ...... 49 3 4 Migration energies and and effective frequencies for the vacancy assisted vanadium self diffusion mechanism for the three ordered V 2 C structures calculated using DFT. ................................ ................................ .... 54 3 5 Correlation factors for vacancy assisted sel f diffusion in an hcp lattice calculated using values of obtained for vanadium self diffusion in the 2 C ................................ ................................ ................................ ...... 55 3 6 Activation energies a nd diffusion prefactors for vacancy assisted vanadium self diffusion for the three ordered V 2 C structures calculated using DFT. ................................ ................................ ................................ ........................ 56 3 7 Migration energies and and effective frequencies for carbon self diffusion for each of the three ordered V 2 C structures calculated using DFT. ................................ ................................ ................................ ........................ 61 3 8 Activation energies and diffusion prefactors for carbon self diffusion for each of the three ordered V 2 C structures calculated using DFT. ......................... 62 4 1 Theoretical equilibrium lattice constants, cohesive energies and magnetic moments for ground state phases of V, Fe, Ni, Ce, Nd and U. ............................ 71 4 2 Supercell dimensions and calculation parameters used in defect calculations for V 2 C ................................ ................................ ................................ ................ 73 4 3 Substitutional and interstitial defect formation energies and impurity volumes. ... 76 4 4 Impurity vacancy binding energies for each of the six in plane and six out of plane nearest neighbor s ites on the vanadium sublattice ................................ .... 79 5 1 Calculated elastic constants and glass transition temperature T g for Cu 30 Zr 70 Cu 50 Zr 50 and Cu 70 Zr 30 metallic glass samples. ................................ ...................... 91 5 2 Calculated parameters for both the elastic and plastic branches of the U p U s Hugoniot for Cu 30 Zr 70 Cu 50 Zr 50 and Cu 70 Zr 30 metallic glass. .............................. 100
10 7 1 Comparison of equilibrium lattice constants for the low temperature phase of LaBGeO 5 ................................ ................................ ................................ .............. 134 7 2 Comparison of atom positions for the low temperature phase of LaBGeO 5 ..... 136 7 3 Comparison of B O, Ge O and La O interatomic distances (in ) for the low temperature phase of LaBGeO 5 ................................ ................................ ......... 137 7 4 Comparison of atom positions for the high temperature phase of LaBGeO 5 .... 141 7 5 Comparison of B O, Ge O and La O interatomic distances (in ) for the high temperature phase of LaBGeO 5 ................................ ................................ ......... 143 7 6 point modulation of the high temperature phase and those between the initial and saddle points along the MEP. ............ 148
11 LIST OF FIGURES Figure page 1 1 Dependence of fracture strength on strain rate for several metall ic glasses ...... 22 3 1 Symmetry unique superstructures generat ed from a 4 site hexagonal cell in order of increasing volume. ................................ ................................ .................... 40 3 2 Energies for a subset of low energy V 2 C structures enumerated from derivative superstructures with volumes times the volume of the base 4 site hexagon al cell. ................................ ................................ ................................ 42 3 3 Unit cells of three ordered V 2 C structures including the low energy orthorhombic structure belonging to space grou p Pnnm ................................ ................................ ................................ ................. 43 3 4 2 C from DFT calculations. ................................ ................................ ................................ ............ 44 3 5 the low energy orthorhombic structure belonging to space group Pnnm the 2 C and the 2 C. ................................ ................................ ........................ 46 3 6 Theoretical phonon DOS for the low energy orthorhombic structure belonging to space group Pnnm 2 2 C. ................... 47 3 7 Calculated vacancy formation energies as a function of the number of atoms in the supercell for the low energy structure belonging to sp ace group Pnnm 2 2 C. ................................ ........................ 48 3 8 I n plane and out of plane jumps for vanadium to a nearest neighbor site within the vanadium sublattice without octahedral sites shown and with octahedral sites shown. ................................ ................................ .......................... 51 3 9 Theoretical MEPs for vanadium vacancy migration for in plane and out of I NEB calculations. ................................ .. 53 3 10 Coordination polyhedra for the octahedral and tetrahedral interstitial sites in V 2 C. ................................ ................................ ................................ ......................... 57 3 11 Octahedral octahedral, labeled A and B, and octahedral tetrahedral octahedral, labeled C, site jumps considered as likely migration paths for carbon self diffusion in V 2 C without and with the vanadium atoms shown. ......... 58 3 12 Theoretical MEPs for carbon i nterstitial migration in plane and out of plane NEB calculations. ................................ ........................... 60
12 3 13 2 C and the low energy orthorhombic structure described by space group Pnnm ........... 64 3 14 Calculated vanadium self diffusion coefficients for the three ordered V 2 C structures and experimental self diffusion coefficients for Nb in NbC x and Ti in TiC x ................................ ................................ ................................ .................... 66 3 15 Calculated carbon self diffusion coefficients for the three ordered V 2 C structures and experimental self diffusion coefficients for C in NbC 0.868 and C in V 6 C 5 ................................ ................................ ................................ .................... 6 7 4 1 Unit cell of low V 2 C. ........................ 72 4 2 Convergence of substitutional and interstitial defect formation energies with respect to the size of the defect supercell ................................ ............................ 75 4 3 Impurity V 2 C structure: six in plane pairs have indices 1 6 and out of plane pairs have indices 7 12. ................................ ..................... 78 4 4 Average in plane and out of plane impurity vacancy binding energies and distance from the ideal vanadium site for substitutional impurities X = Fe, Ni, V 2 C. ................................ ................................ .......................... 80 4 5 MEP s for migration of Fe and Ni impurities for each of the six in plane and six out of plane vanadium sites. ................................ ................................ ............ 82 5 1 Zero temperature stress strain curves and sound speed curves calculated using the Cheng potential for Cu Zr( Al) and the Mendelev potential for Cu Zr. ................................ ................................ ................................ ............................ 87 5 2 MSDs for liquid and metallic glass and Cu Cu RDFs for liquid and metallic glass. ................................ ................................ ................................ ....................... 90 5 3 Volume temperature relationships during cooling for Cu 30 Zr 70 Cu 50 Zr 50 and Cu 70 Zr 30 systems. ................................ ................................ ................................ ... 91 5 4 Shock wave profiles of shock pressure and shear stress from MD simulations of Cu 50 Zr 50 metallic glass at U p = 0.125, 0.5 and 1.5 km/s. ................................ .. 94 5 5 Shock pressure at the HEL P HEL as a function of final shock pressure for Cu 30 Zr 70 Cu 50 Zr 50 and Cu 70 Zr 30 metallic glass. ................................ ..................... 96 5 6 Shock pressure distributions at different times and positions for Cu 50 Zr 50 metallic glass at U p = 1.5 km/s and U p = 0.5 km/s. ................................ ............... 97 5 7 Shock velocity U s as a function of piston velocity U p for Cu 30 Zr 70 Cu 50 Zr 50 and Cu 70 Zr 30 metallic glass. ................................ ................................ ................... 99
13 5 8 Shock wave profiles at t = 100 ps for Cu 30 Zr 70 Cu 50 Zr 50 and Cu 70 Zr 30 metallic glass samples at U p = 1.5 km/s. ................................ ................................ ........... 101 5 9 Shock induced temperatures as a function of shock pressure. .......................... 102 5 10 Peak shear stress and flow stress as a function of the final shock pressure. .... 103 5 11 Cu Cu RDFs in the aftershock flow of Cu 50 Zr 50 metallic glass for four different shock intensities. ................................ ................................ ................................ .. 105 5 12 Shock pressure at the HEL P HEL versus mean pressure P m at the HEL. ........... 106 6 1 Pressure P xx and shear stress profiles at t = 125 ps for Cu 50 Zr 50 metallic glass at different piston velocities. ................................ ................................ ....... 115 6 2 Snapshots of von Mises shear strain for Cu 50 Zr 50 metallic glass at different piston velocities ................................ ................................ ................................ ... 116 6 3 Snapshot at t = 40 ps of von Mises shear strain i vM and profiles of vM (x) and shear stress for Cu 50 Zr 50 and U p = 0.375 km/s. ................................ ................... 118 6 4 Peak shear stress peak as a function of the mean pressure P m and flow stress flow as a function of the final shock pressure P xx ................................ ..... 119 6 5 Snapshots of von Mises shear strain for Cu 30 Zr 70 Cu 50 Zr 50 and Cu 70 Zr 30 metallic glasses at P xx = 20 GPa. ................................ ................................ ........ 121 6 6 Average von Mises shear strain in the quasi steady flow state B behind the shock front as a function of the final shock pressure P xx ................................ ... 122 6 7 Diffusivity as a function of shock pressure P xx for Cu Zr metallic glasses from this work and Cu 46 Zr 54 metallic glass from Hugoniostat simulations. ................. 123 6 8 Representative Cu centered clusters corresponding to t he given set of Voronoi indices ................................ ................................ ................................ .... 125 6 9 Fraction of Cu centered Voronoi polyhedra at the quasi steady flow state behind the shock front as a function of the pisto n velocity U p ........................... 126 6 10 Fraction of Cu centered Voronoi polyhedra that have been transformed into other types after transitioning from the unshocked state to the quasi steady flow state B behind the shock front ................................ ................................ ..... 128 7 1 Structure of the low temperature phase of LaBGeO 5 from the side and from the top. ................................ ................................ ................................ .................. 133 7 2 Structure of the LaO 9 polyhedra with attached GeO 4 (green) and BO 4 (blue) tetrahedra and without. ................................ ................................ ......................... 135
14 7 3 Electronic band structure and DOS for the low temperature phase of LaBGeO 5 calculated usin g DFT. ................................ ................................ .......... 138 7 4 Phonon dispersion and DOS for the low temperature phase of LaBGeO 5 calculated using DFT. ................................ ................................ ........................... 139 7 5 Structure of the high temperature phase of LaBGeO 5 ................................ ....... 140 7 6 Structure of the LaO 9 polyhedron from the low temperature phase of LaBGeO 5 and the LaO 10 polyhedron from the high temperature phase. ........... 142 7 7 Phonon dispersion and DOS for the high temperature phase of LaBGeO 5 calculated using DFT. ................................ ................................ ........................... 144 7 8 Energy as a function of the relative displacement of atoms from the high temperature phase of LaBGeO 5 M and K points. ................................ ................................ ................................ ..... 145 7 9 MEP between the low temperature (ferroelectric) and high temperature (paraelectric) phases of LaBGeO 5 and magnitude of polarization along the c axis at different points along the MEP ................................ ................................ 146 7 10 Illustration of the rigid rotation of the BO 4 tetrahedra occurring during the ferroelectric phase transition in LaBGeO 5 ................................ .......................... 149
15 LIST OF ABBREVIATIONS CI NEB Climbing image nudged elastic band DFT Density functional theory DOS Density of states EAM Embedded atom method FCCI Fuel cladding chemical interaction GGA Generalized gradient approximation HEL Hugoniot elastic limit HNF Hermite normal form LAMMPS Large scale Atomic/Molecular Massively Parallel Simulator LDA Local density approximation MEP Minimum energy path MG Metallic glass MSD Mean square displacement NN Nearest neighbor PAW Plane augmented wave PBE Perdew Burke Ernzerhof PES Potential energy surface RDF Radial distribution function SFR Sodium fast reactor SRO Short range order STZ Shear transformation zone TST Transition state theory UFF Universal force field VASP Vienna Ab Initio Simulation Package
16 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 TRANSPORT AND PHASE TRANSITIONS IN CERAMICS AND AMORPHOUS METALS BY ATOMIC LEVEL SIMULATION By Brian Demaske December 2017 Chair: Simon R. Phillpot Major: Materials Science and Engineering Microscopic phenomena of diffusion, phase transitions and high strain rate mechanical behavior are investigated in ceramics and amorphous metals using atomic scale simulation methods. The defect and diffusion behavior of V 2 C is investigated using density f unctional theory calculations owing to its potential use as a diffusion barrier in nuclear fuel cladding applications Self diffusion coefficients for carbon and vanadium are calculated using methods based on transition state theory and compared to existin g experimental measurements Dilute m etal impurities are found to p refer substitutional vanadium sites over interstitial sites Impurity vacancy binding is found to increase with impurity size with the largest impurities forming clusters with neighb oring v anadium vacancies Shock compression of Cu x Zr 100 x (x = 30, 50 and 70) metallic glasses is carried out over a wide range of shock intensities using molecular dynamics simulations An elastic plastic response characterized by a narrow elastic precursor and sluggish plastic shock wave is observed for all compositions. The Hugoniot elastic limit of the simulated metallic glasses is found to be strongly pressure dependent, in contrast to crystalline metals. Yielding of the metallic glass occurs via nucleation of shear transformation
17 zones with a gradual shift towards homogeneous flow at high shock intensities. Local s hear resistance of the metallic glass is correlated with the degree of icosahedral order of Cu centered clusters which is f ound to increase with Cu content. The shock ind uced melting transition is found to occur at higher shock pressures for samples with higher Cu content Ferroelectricity in LaBGeO 5 with stillwellite structure is investiga ted using density functional theory calculations. Phonon dispersion calculations of the high temperature paraelectric phase reveal an unstable mode that is shown to drive the ferroelectric phase transition. The main feature of the ferroelectric distortion is a rigid rotation of BO 4 t etrahedra which agrees with previous short range potential calculations A spontaneous polarization of 4.6 C/cm 2 is calculated for LaBGeO 5 using the modern theory of polarization. This value lies within the experimental range of 2.7 12 C/cm 2
18 CHAPTER 1 INTRODUCTION 1.1 General Overview On the smallest scales, the behavior of a material is ultimately determined from the interactions among its constitu ent atoms. The potential energy of a system comprised of N atoms is a 3 N dimensional function of the atomic coordinates. All possible sets of atomic coordinates, or simply configurations, make up the potential energy surface (PES) of a given system. Along the PES, there exist local minima that correspond to configurations of stable or quasi stable structures, such as crystal phases, glasses or molecules. A transition from one distinct local potential energy minima to another occurs along a minimum energy pa th (MEP) of the PES that passes through a saddle point or transition state. At the transition state, the potential energy attains a local minimum in all but one dimension in which it attains a local maximum. In addition to the potential energy, atoms in re al materials also possess kinetic energy due primarily to thermal motion. At low temperatures, systems can only sample a small area of the PES around a local energy minimum. As the temperature is increased, however, the number of configurations that the sy stem passes through increases, which gives it the ability to sample more of the PES. Occasionally, the system will overcome an energy barrier between two nearby energy minima by passing through a transition state. Together the height of the energy barrier and the temperature of the system dictate the rate of the process: lower barriers and higher temperatures correspond to faster processes. Many phenomena, such as phase transitions, diffusion, chemical reactions, glass formation, structure prediction, etc., can be easily understood on the basis of the PES or
19 energy landscape paradigm. For example, the process of glass formation is the result of cooling a liquid at a rate that is faster than the competing crystallization rate  While traversing the PES, the system becomes trapped in a lo cal energy minimum corresponding to an amorphous or glassy configuration and due to insufficient thermal energy is unable to overcome the barrier to the lower energy crystalline state. The idea of the energy landscape has proven extremely useful in atomic scale simulations in which the configurations of and interactions among individual atoms are controlled and specified by the user. In this work, we use both static and dynamic simulation methods to study problems related to lattice diffusion, phase transi tions and plastic deformation processes in ceramics and metals. Static simulation methods are used to search out energy minima on the PES to determine stable ordered crystalline structures of V 2 C and transition state searches are conducted between neighbor ing energy minima to determine MEPs for diffusion. The dynamics of the ferroelectric phase transition in the stillwellite compound LaBGeO 5 are determined by characterizing the PES around a known transition state configuration Dynamic methods are used to s imulate the trajectories of atoms within a metallic glass (MG) under the action of high strain rate shock compression. MGs correspond to quasi stable energy minima along the PES P lastic deformation in MGs has been shown to be a thermally activated process whereby a small cluster of atoms undergoes a collective shear mo tion in response to an applied shear stress  1.2 Motivation of Problems 1.2.1 Vanadium Carbide Diffusion Barrier Nuclear fuel rods consist of fuel pellets encased within a metallic tube or cladding. During operation, the gap between the fuel and cladding closes and
20 components of the fuel, including fission products, and cladding can diffuse into one another forming an interaction layer [3 6] This process is known as f uel cladding chemical interaction (FCCI) and can result in the formation of brittle or low melting point phases that thin the cladding wall and make it susceptible to failure. Next generation nuclear reactors, such as sodium cooled fast reactors (SFRs) are being developed that enable higher fuel burn ups and better thermal efficiency when compared to current generation reactors. In SFR designs, U Zr or U Pu Zr metallic alloys represent typical fuels and ferritic/martensitic steels represent typic al cladding materials  The metal metal interface can exacerbate FCCI in contrast to the metal oxide interface present in previous generation reactors. Moreover the high service temperatures and irradiation levels in SFRs can also act to further enhance FCCI [8 11] so new methods are being sought to reduce the interdiffusion of fuel and cladding components. Recently, va nadium has emerged as a promising candidate material to act as a diffusion barrier between the fuel and cladding [12 14] To en sure there is good bonding between the diffusion barrier and cladding, a pack cementation diffusion coating (PCDC) process was developed to deposit vanadium on HT 9 steel  a potential candidate cladding material in design of SFRs  Decarburization of the steel during the PCDC process leads to the formation of a thin layer with nearly V 2 C stoi chiometry. O ut of pile d iffusion couple experiments between Ce acting as a fuel surrogate, and HT 9 steel with and without the V 2 C coating layer were performed and i t was found that the interdiffusion layer is effectively eliminated when the vanadium coa ting i s present 
21 The kinetics of diffusion in V 2 C is largely unknown M ost of the available experimental diffusion data has been obtained for carbon self diffusion in the cubic transition metal carbides [16 26] Self diffusion data on metal atom self diffusion is rare [27,28] and experiments on impurity diffusion are nonexistent. Moreover, little theoretical wo rk has been done to address lattice diffusion in the transition metal carbides. The goal of Chapter s 3 and 4 of this work is to characterize the diffusion behavior of intrinsic defects and impurities in V 2 C. Self diffusion of vanadium and carbon will be treated in Chapter 3 guided by available experimental results for self diffusion in the cubic transiti o n metal carbides. Chapter 4 will proceed to address the behavior of impurities in V 2 C specifically addressing those impurities originating in nuclear fuel and steel cladding materials Calculations are performed using density functional theory ( DFT ) as it has proven effectiveness in accurately describing self diffusion and impurity diffusion in metals [29 38] and can treat a wide variety of impurity types at no added cost. 1.2.2 High strain rate Deformation of Metallic Glasses In terms of mechanical properties, MGs outperform crystalline alloys in yield strength, elastic strain limits and energy storage [2,39 49] Their superior performance is grounded in t heir amorphous structure, which leads to a fundamentally different plastic deformation mechanism than the dislocation mediated plastic response of metallic crystals. Quasistatic compression and tension experiments at room temperature have shown that plastic deformation in MGs occurs inhomogeneously through the nucleation and growth of shear bands [44,2,46 48,50] Under tension, only a single dominant shear band may be activated resulting in catastrophic failure of the MG with little to no macroscopic plasticity. However, under constrained loading
22 geometries, like uniaxial compression or bending, serrated flow i s ob served and some macroscopic plasticity (~1% plastic strain) is achieved before failure The serrated flow correspond s to emission of localized shear bands in a step by step manner. Figure 1 1. Dependence of fracture strength on strain rate for several m etallic glasse s. Reprinted from  with permission from Elsevier. Less is known about the deformation behavior of MGs under d ynamic loading conditions The f racture strength s of several representative MGs under uniaxial stress conditions are shown in Fig. 1 1 for a range of strain rates from 10 4 to 10 4 s 1 The b ulk (> 1 mm thickness) MGs were found to exhibit rate independent behavior or a sligh tly negative response with increasing strain rate At the highest strain rates, the fracture surfaces of the bulk MG samples are rough due to simultaneous ope ration of multiple shear bands [2,47,51,52] The large drop in strength for the Zr 65 Al 10 Ni 10 Cu 15 MG ribbon (0.05 mm) is presumed to be due to surface effects. Although there is some agreement
23 among the different studies the effect of strain rate on the yield behavior of MGs has not been shown d efinitively Shock waves have long been used as a tool to measure the strength of materials at very high strain rates (> 10 5 s 1 ) [53 55] In ex periment s high velocity projectiles or laser pulse s are used to generate shock waves which propagate through the target at supersonic velocities. Material entering the shock is rapidly compressed along the shock wave propagation direction Such compression results in a state of uniaxial strain within the sample, where as the results discussed above were obtained for uniaxial stress states. If the amplitude of the shock wave exceeds the Hugoniot elastic limit (HEL), the n the material begin s to yield. Shock compression experiments have been performed on several MGs [56 66] However, d ue to the high speed of the shock wave experimentalists are limited in their ability to obtain thermodynamic or microstructural information in sit u Time resolved X ray diffraction measurements have been useful in analyzing shock deformation in crystals [67 71] but have yet to be applied to MGs. Molecular dynamics (MD) simulations have been used to study shock deformation in many different types of materials [72 93] Results of these simulations have provided insights into the microscopic defo rmation behavior of single and nanocrystalline metals [75,83,87 89,94] as well as resolved fine details of shock wave structure [78,85,86] In contr ast to experiments, MD simulations provide the complete atomic scale dynamics of the system starting from the initial loading event to reflection of the shock at the opposite surface Chapters 5 and 6 of this work will characterize the shock response of Cu Zr MGs using MD simulations Specific focus will be placed on the thermodynamic response of the MG as well as changes to the microstructure that
24 occur during shock compression. Such simulations will provide insights into the behavior of MGs under shock lo ading as well as the microscopic plastic deformation mechanism at high strain rates. 1.2.3 Lanthanum Borogermanate G lass ceramics The l anthanum borogermanate La 2 O 3 B 2 O 3 GeO 2 system has good glass forming ability  Wh en subjected to heat treatment, some glass compositions u ndergo monophase crystallization into LaBGeO 5 [95 97] which is isomorphous with stillwellite CeBSiO 5 At room temperature, LaBGeO 5 is ferroelectric [98 102] and, when doped with Nd 3+ or Pr 3+ LaBGeO 5 can be used a self frequency doubling laser [103,104] Surface and bulk crystallization of lanthanum borogermanate glasses can be achieved through various heat treatment processes. By adjusting the temperature and duration of the heat treatment, it is possible to genera te transparent glass ceramic composites with nanometer sized crystallites of LaBGeO 5 [105 109] Such transparent glass ceramic nanoc omposites have applications in optoelectronics, integrated optics, photonics and communication technologies due to the nonlinear optical properties of the embedded LaBGeO 5 crystals. In addition to the heat treatment method discussed above it has been shown that crystallization of lanthanum borogermanate glasses is also possible through focused laser irradiation [110 115] Heating of the glass by a pulsed femtosecond laser leads to localized melting of the glass and formation of a LaBGeO 5 seed crystal at the melt glass interface. By continuously moving the s ample during laser irradiation crystalline tracks can be grown within the bulk or on the surface of the glass Moreover, provided the bending angles along a track are gradual, the ori entation of the crystal varies continuously through the bend without forming grain boundaries 
25 Interestingly, the polar axis of the LaBGeO 5 crystal aligns with the direction of the laser writing. Recently, it has been shown that these crystalline tracks can function as waveguide s with second harmonic generation properties  The Curie temperature for LaBGeO 5 is ~530 C [98,99] Below this temperature, LaBGeO 5 attains a nonzero spontaneous polarization along the c axis that can be reversed u nder an electric field. Measurements of the spontaneous polarization at room temperature have a range of 2.7 12 C/cm 2 [98 102] Neutron diffr action data show continuous variation of the cell volume through the temperature range 200 700C  which is evidence for the second order na ture of the phase transition at 530C Subsequent experimental measurements  and force field modelling  of the vibrational properties of LaBGeO 5 found the ferroelectric phase transition is dominated by rotation of BO 4 tetrahedra. However, the lattice dynamics of the ferroelectric phase transition could not be conclusively shown from these works Chapter 7 o f th is work will use accurate DFT calculations to characterize the low and high temperature phases of LaBGeO 5 and to determine the complete lattice dynamics of the ferroelectric phase transition From these results, the spontaneous polarizatio n of ferroel ectric phase is calculated within the framework of the modern theory of polarization [1 19 123] 1.3 Outline This dissertation is organized as follows. In Chapter 2, brief o verviews of DFT and MD are given emphasizing their specific applications to the calculations performed in this work. In Chapter s 3 and 4, the DFT calculations of the energetics and diffusion behavior of intrinsic defects and select impurities in V 2 C are presented The focus of Chapter 3 is on the intrinsic defect and self diffusion behavior. Self diffusion coefficients are calculated under the assumptions of a vacancy assisted mechanism for vanadium
26 and an interstitial mechanism for carbon. In Chapter 4 the site preference and impurity vacancy binding energies are determined for select impurity types as a first step towards characterizing the impurity diffusion behavio r. Impurity types considered are Fe and Ni, U and Ce and Nd, which represent major components of the clad, nuclear fuel and fission products respectively Migration barriers are calculated for the two smallest impurities, Fe a nd Ni. In Chapters 5 and 6, M D simulations of shock compression of Cu Zr MGs are presented The focus of Chapter 5 is on the shock wave structure, Hugoniot states and yield criterion. In Chapter 6, the microstructural changes of the MGs during shock loading are investigated with empha sis on the influence of short range order (SRO) on the plastic deformation behavior. In both Chapters 5 and 6, the composition of the MG is varied to determine its effect on the shock response In Chapter 7 DFT calculations o f the ferroelectric phase tran sition and spontaneous polarization in LaBGeO 5 are presented. A collected summary of results and ideas for future work are given in Chapter 8.
27 CHAPTER 2 OVERVIEW OF SIMULATION METHODS D efect and diffusion behavior in V 2 C as well as the ferroelectric phase transition and spontaneous polarization in LaBGeO 5 will be investigated using first principles DFT calculations. Shock compression of Cu x Zr 100 x MGs will be simulated using MD 2.1 Density Functional Theory DFT can treat upwards of a thousand atoms with routin e calculations being performed on systems with several tens to a few hundred atoms. Such small system sizes preclude direct calculation of diffusion coefficients in solids, because of the long time and finite temperatures needed to sample diffusive process es. Fortunately, diffusion problems can be readily addressed via static DFT calculations using ideas from transition state theory (TST)  To use TST methods, the diffusion mechanism mu st be known ahead of time, i.e., the initial and final states of the diffusing atom must be known. As such, only simple diffusion problems can be addressed in this way. A method for calculation of diffusion coefficients in simple crystal structures using D FT is given in Ref.  Macroscopic polarization can also be treated within the framework of DFT owing to the development of the modern theory of polarization [119 123] D FT is a reformulation of nonrelativistic quantum mechanics in the language of density functio nals The f oundation of DFT is the two Hohenberg Kohn theorems  The first theorem states that the external potential of an N electr on system is a unique functional of the ground state electron density By implication, a ll other observa bles, including the total energy, are functionals of The s econd theorem states that where is some other density satisfying t he nec essary boundary conditions Thus, the ground state energy can be obtained by minimizing
28 K ohn and Sham showed that the minimization problem can be mapped to a set of single particle Schr dinger like equations, known as the Kohn Sham equations  w here are the single particle orbitals of the non interacting system. The Hartree potential represents the Coulomb repulsion between particles All non classical interactions are contained in the exchange correlation potential The e xact form of the exchange correlation functional is unknown, so approximations must be made T he two most common approximations used in the study of condensed matter systems are the local density approximation (LDA) and the generalized gradient approximation (GGA). In the LDA, where is the exchange correlation energy of a homogeneous electron gas o f density The GGA expands on the LDA by including the g radient of the electron density within the calculation of There are many different variants of GGA depending on the form of In addition to the Perdew Zunger parameterization of the LDA  this work u ses the Perdew Burke Ernzerhof (PBE) form of the GGA  and a variant designed specifically for solids, PBEsol  With an appropriate expression for t he Kohn Sham equations can then be solved self consistently to determine and, hence, the ground state energy Note that given the exact form of the charge density is equal to for the true N electron system. In practice, are expressed in terms of either localized (atomic like) basis sets, like Gaussians, or delocalized basis sets, like plane waves. DFT c od es, such ( 2 1 )
29 as VASP  that use p lane wave basis sets are partic ularly well suited to solve problems in periodic systems, which is one of the main reasons for their popularity among researchers in condensed matter physics and materials science In DFT applications is specified fro m the positions of atoms or nuclei, which are assumed to remain fixed when solving for the electronic ground state In a nutshell, this is the Born Oppenheimer or adiabatic approximation. Once the electronic ground state is obtained self consistently, i nteratomic forces and stres ses can be obtained via t he Hellman Feynman theorem By minimizing t he interatomic forces and stresses in the cell one can calculate the equilibrium structure of the system. In addition, by considering finite changes in atom positions and/ or cell vectors, one can also calculate elastic constants, phonons and other aspects of the PES as well as a range of other properties 2.2 Modern Theory of Polarization With the advent of the modern theory of polarization [119 123] calculation of the macroscopic polarization is now a standard feature of most DFT codes. According to the theory, the macrosc multivalued quantity defined modulo a vector known as the polarization quantum, where is the electronic charge and ( i = 1, 2, 3) is a lattice vector. As a consequence, the polarization for a given system is not uniquely defined, but is instead represented by a lattice of values with spacing given by  Moreover, even nonpolar systems can exhibit nonzero polarization equal to or where is an integer  Different values of correspond to different branche s within the polarization lattice. The spontaneous polarization of a ferroelectric is defined as the
30 polarization difference between a polar system and a reference nonpolar system or, equivalently, as one half the polarization difference between two polar systems with enantiomorphous structures  In practical calculations of the spo ntaneous polarization, it is useful to calculate the polarization at several points along a continuous transformation path between these two systems to ensure that the two end structures reside on the same branch of the polarization lattice. 2.3 Molecular Dyna mics Material entering a shock front is compressed uniaxially at very high strain rates, which leads to heating and stress generation. For moderate shocks, i.e., those that occur in most materials applications, the temperatures and stresses are low enough to safely ignore electronic effects. Because of the high strain rates involved in shock compression of solids, it is necessary to explicitly treat interactions among individual atoms to understand how the material will respond near the shock front Thus, M D has become the natural choice for simulating shock waves in solids [72 93] In MD simulation the classical equations of motion are solved numerica lly to yield the trajectories for a system of N interacting particles, usually atoms Given a set of boundary and initial condit ions, positions and velocities particle traject ories are updated after a time step using a finite difference algorithm. There are many choices of algorithms, Verlet  velocity Verlet  predictor corrector  etc. In this work, the velocity Verlet algorithm is used, in which the position and velocit y of a particle at time are given by ( 2 2 )
31 where is the acceleration. From t he acceleration of each particle is proportional to the net force ac ting on it where is the potential energy function or interatomic potential, of the system and After the positions and velocities of all particles have been updated using Eq. 2 2, the simulation time is advanced and the process is repeated for the next time step. During a n MD run pr operties for individual parti cles can be calc ulated using relations from statistical mechanics For example, the instantaneous temperature for a particle with three degrees of freedom By perf orming averages over a group of particles or a series of time steps or both, one can extract thermodynamic properties of the system, such as kinetic and potential energies, temperature, pressure tensor, flow velocity, density etc. Other information about the system can be determined by monitoring t ra jectories, like the mean square displacement (MSD) or atomic strain tensor, or by looking at the local environments around particles, like radial distribution functions (RDFs ) or Voronoi polyhedra. Proper solutions to equations of motion wi ll conserve the total energy ( E ), volume ( V ) and number of particles ( N ) in the system. Such simulations correspond to the NVE or microcanonical ensemble. It is possible to modify the equations of motion to simulate other thermodynamic ensembles. Some of t he most common ensembles used in MD simulations other than NVE are the isochoric canonical ensemble ( NVT ) and the isobaric canonical ensemble ( NPT ). In this work, both the NVE and NPT ensembles are used. In the NPT ensemble, the temperature and pressure of the system are controlled by a thermos tat and barostat, respectively Th e
32 NPT calculations in this work are performed using a Nos Hoover thermostat [136 138] and Par r inello Rahman barostat  MD simulations presented in this work are performed using LAMMPS  The interatomic potential determines the forces on particles and, therefore, the dynamics of the system under st udy. In ab initio MD, the potential is supplied directly from electronic structure calculations, like DFT. However, after each time step, the potential must be recalculated as the atoms move and the external potential changes. Therefore, ab initio MD is computationally demanding and is limited to small systems (~ 1 1 0 00 atom s) and short time durations (~1 100 ps). In contrast, in classical MD electron ic degrees of freedom are not treated explicitly and the interatomic potential is fixed ahead of tim e, which allows for much larger system sizes and longer simulation times. The disadvantage is that the potential must be carefully chosen for the problem at hand with no guarantee of transferability. Many types of potentials exist for many different types of systems that range in complexity from simple pair potentials, like Lennard Jones  to complex reactive force fields, like ReaxFF  2.4 Embedded Atom Method T he MD simulations in this work are concerned with metallic alloys, specifically Cu Zr MGs. P erhaps the most widely used and successful approaches to modeling interatomic potentials in metals is the embedded atom method (EAM) [143,144] The potential energy in the EAM formalism is given by ( 2 3 )
33 where is a pair potential between atom types and and is an embedding energy for atom The charge density at site is where is the charge density associated with atom In the original derivation, is a measure of the local electron density of each atom and represents the energy associa ted with placing an atom in a region with electron density The pair potential represents short range interactions between cores. However, most current implementations treat each of the terms in Eq. 2 3 as fitting functions and the potential is fit to a database comp ri sing experimental and ab initio data  The success of EAM potentials lie in their ability to accurately model metallic bonding through the many body embedding term
34 CHAPTER 3 INTRINSIC DEFECTS AND SELF DIFFUSION IN VANADIUM CARBIDE 1 3.1 Background Transition metal carbides exhibit a combination of excellent mechanical and thermal properties, such as high hardness, high melting temperature, good wear resistance and high thermal conductivity [146 148] These properties combined with the success of vanadium metal as a cladding l iner for nuclear fue l rods [12,149] have sparked interest in the use of vanadium carbide as a diffusion barrier to help mit igate FCCI in next generation nuclear reactors  Using Ce as a stand in for the metallic fuel and fission products, diffusion couple experiments were performed at 660C (933 K) for 100 hours on HT9 steel with and without the vanadium carbide coating  Results indicate an almost complete elimination of the interdiffusion layer when the vanadium carbide coating is present. The coating process used in Ref.  results in the formation of a vanadium carbide layer having clo se to V 2 C stoichiometry. According to existing phase diagrams [150 152] a disordered V 2 C y phase is stable at high temperatur es and is characterized by a hexagonal vanadium sublattice with carbon atoms randomly occupying y /2 of the available octahedral sites. Ordering of the carbon atoms occurs at lower temperatures resulting in the formation of two ordered phases, V 2 C and V 2 C. Of these two phases, V 2 C is the experimental ground state structure The order disorder transition temperature between V 2 C and V 2 C y is estimated to be in the range of 1100 1600 K  though the exact temperature is unknown. 1 The work described in this chapter has been published in B.J. Demaske, A. Chernatynskiy, and S.R. Phillpot, J. Phys: Condens. Matt. 29, 245403 (2017).
35 S elf diffusion of carbon in cubic transition metal carbides with NaCl type structure has been inv estigated experimentally for NbC  ZrC [17,18] TiC [16,17,22] and V 6 C 5  whereas metal atom self diffusion has been characterized only in NbC x  and TiC x  Experimental diffusion studies on transition metal c arbide s with non NaCl type structures are scarce with only a few studies on WC  Mo 2 C  and Fe 3 C  Furthermore, experiments have yet to address the diffusion of impurities in any of the transition metal carbides, including V 2 C. From a theoretical perspective, only basic studies of the structural, electronic and elastic pr operties of V 2 C ha ve been performed [155 157] No work has yet been done to address the diffusion of intrinsic defects or impurities in V 2 C. The objective of the following two chap ters will be to investigate intr insic defects and impurities in V 2 C using first principles DFT calculations This chapter will address the intrinsic defect behavior and self diffusion of carbon and vanadium while Chapter 4 will address the energetics and kinetics of dilute metal impurities representative of the fuel clad system Self diffusion experiments on the transition metal carbides found that metal atom diffusivity is independent of carbon concentration, whereas carbon diffusivity is not. Therefore, i t can b e inferred from these results that carbon and metal atom diffusion proces ses are uncorrelated. R esearchers have postulated the most likely self diffusion mechanisms in NaCl type M x C y including (1) the direct migration of an atom to a nearest neighbor (NN) vacant site and (2) a two jump sequence in which the migrating atom temporarily occupies a tetrahedral site before advancing to a NN vacant site. The first mechanism is thought to occur primarily for metal atoms [27,28] while the second
36 mechanism is thought to occur for carbon [20,24] The self diffusion behavior of binary transition metal carbides has been summarized in Refs.  and  This chapter will begin to fill the gap in understanding self diffusion in non cubic transition metal carbides from a theoretical perspective using V 2 C as a case study owing to its potential use as a cladding liner in nuclear fuel rods Expected a pplication temperatures are low enough such that carbon or dering is likely. Therefore, the focus is o n ordered V 2 C structures. An exhaustive search of different possible ordered structures up to a certain specified volume and satisfying the hexagonal symmetry of the vanadium sublattice is conducted to determine the lowest energy configurations maintaining V 2 C stoichiometry. Equilibri um structural, electro nic and elastic properties of representative V 2 C structures are calculated and compared to existing experimental and theoretical data prior to intrinsic defect and self diffusion calculations. Predicted carbon self diffusion coefficie nts, diffusion prefactors and activation energies are com pared to available experiments for the disordered cubic phases of V 6 C 5  and NbC 0.868  while theoretical results for vanadium self diffusion are compared to available experiments f or self diffusion of Nb in NbC x  and Ti in TiC x  Experimental diffusion data is measu red for single crystals over the temperature range 2000 2500 K. 3.2 Simulation Methods The equilibrium structural, electronic and elastic properties and intrinsic defect energetics and self diffusion coefficients of V 2 C a re calculated using VASP [131,159] which employs the plane augmented wave ( PAW ) method to simplify treatment of core electrons [160,161] In our calculations, 4 electrons for C (2 s 2 2 p 2 ) and 13 electrons for V (3 p 6 3 d 3 4 s 2 ) a re treated as valence. For comparison purposes, both the PBE form of the GGA  and the Perdew Zunger parameterization of t he LDA  a re used to
37 approximate the exchange centered Monkhorst Pack k point meshes for each supercell a re generated using the fully automatic scheme implemented in VASP in which the number of subdivisions in the k point mesh along a reciprocal lattice vector is determined according to ( 3 1 ) where l is a parameter that can be varied to change the density of the k point mesh. A plane wave energy cutoff of 600 eV and k point mesh with l = 30 a re found to give energies converged to within 5 meV/atom. Very fine k point meshes with l = 60 a re used in calculations of the el ectronic density of states ( DOS ) and strings of k points along major dir ections in the Brillouin zone a re used in ca lculations of the electronic band structur e. Wave function optimization is truncated when the energy difference between successive electronic steps drops below 10 6 eV. Partial orbital o ccupancies a re treated using the method of Methfessel Paxton (1 st order) with a smearing width of 0.2 eV. Structure relaxation calculations a re carried out using the conjugate gradient algorithm  until t h e maximum force on each atom i s < 0.01 eV/ and th e maximum stress component on the cell i s < 0.01 eV/ The ordered structures, described in the first section of the Simulation Results have different cell shapes and volumes. To test the convergence of the total energies of the ordered structures with respect to number of k points, the k point mesh i s systematically increased from 999 to 1 71717, where a uniform mesh i s used instead of the automatic sche me outlined above. The use of the uniform mesh is preferred as the lattice vectors and, hence, reciprocal lattice vectors vary from structure to structure which would lead to different k point meshes for different structures using
38 the automatic generation scheme The ordered structures a re first relaxed to equilibrium and then a dditional static calculat ions a re performed using the linear tetrahedron method with Blchl corrections  to obtain total energies. Upon increasing the size of the k point mesh from 131313 to 151515, the maximum change in the relative energi es for all ordered structures drops below 0.25 meV. Equilibrium elastic constants a re calculated by considering both finite distortions of the cell with rigid ions  and contributions from ion relaxation  The plane wave energy cutoff and size of the k point mesh a re systematically increased to test convergence of the elas tic constan ts, where a uniform k point mesh i s used instead of the automatic generation scheme as VASP changes the k point set on the fly when deforming the cell Use of the automatic scheme may cause changes in the k point mesh due to the finite deformations of the cell lattice vectors, which can lead to inaccuracies i n the calculated elastic constants Above a plane wave energy cutoff of 1050 eV and a 111111 k point mesh individual elastic constants a re found to vary by < 0.5 GPa. Analysis of the elastic constan t matrices i s performed using the ELATE software package  Dynamical matrices a re obtained using the finite displacement method implemented in the software package Phonopy  with VASP as the force calculator. Migration barriers for self diffusion a re calculated using the cli mbing image nudged elastic band ( CI NEB ) method [168,169] with a force tolerance of 0.01 eV/. By varying the number of intermediate images used in t he CI NEB calculations, we f i nd that the transition states d o not change appreciably beyo nd five images for vanadium atom migration and seven images for carbon atom migration
39 3.3 Simulation Results 3.3.1 Ordered Structures of V 2 C According to X ray diffraction stud ies of nearly stoichiometric V 2 C [150,152,170] at low temperatures V 2 C forms phase (space group ) stable under carbon group Pbcn ) stable below 1100 phase, the vanadium atoms are hexagonal lattice. The two phases differ predominantly in the ordering of carbon atoms. Calculated equilibrium properties for the two experimentally observed ordered phases are given in Table 3 1 phase in our DFT calculations, confirming experimental results. To see how the arrangement of carbon atoms influences the total en ergy, we perform an exhaustive search over all possible configurations of superstructures of volume n that can be generated from the base 4 site hexagonal cell (space group P 6 3 / mmc ), see Fig. 3 1 where n is a multiple of the base cell volume. In Fig. 3 1, l arge (red) spheres represent vanadium atoms and small (gray) spheres represent octahedral sites. Superstructures a re generated using the Hermite normal form (HNF) routines found in the AFLOW software package  Details of superstructure generation using HNF matrices can be found in Ref.  The symmetry unique superstructures of the base cell are shown in Fig. 3 1 in order of increasing volume. For each superstructure, all configurations satisfying the V 2 C stoichiometry a re enumerated by placi ng carbons at half of the octahedral sites, e.g. a superstructure containing m octahedral sites has configurations. The energy of each configuration i s then calculated using the
40 universal force field (UFF) method  The process i s repeated for n from 1 to 8, phases of V 2 C were among the structures generated for n = 4(8) and n = 3(6), respectively. Duplicate configurations with the same UFF energy a re removed and only 110 configurations with the lowest energies a re retained for structural relaxation using VASP. F igure 3 1. Symmetry unique superstructures generated from a 4 site hexagonal cell (space group P 6 3 / mmc ) in order of increasing volume.
41 Table 3 1. phases of V 2 C and the low energy structure belonging to space group Pnnm (see Fig. 3 1). Space group (phase) Pn nm Pbcn P 3 1m ( ) Lattice type orthorhombic o rthorhombic hexagonal Prototype anti CaCl 2 Fe 2 N Fe 2 N a () 4.534 4.549 4.540 4.495 4.551 4.570 5.008 5.005 b () 5.039 5.735 5.726 5.628 5.735 5.742 c () 2.864 5.031 5.031 4.929 5.032 5.026 4.536 4.551 E 0 (eV/atom) 9.5565 9.5554 9.5496 E ZP (meV/atom) 72.9 72.6 72.6 Reference This work (VASP PBE) This work (VASP PBE) VASP PBE Ref.  CASTEP Ref.  CASTEP Ref.  Expt. Ref.  This work (VASP PBE) Expt. Ref.  Note: theoretical resu lts obtained at zero temperatur e
42 Figure 3 2. Energies for a subset of low energy V 2 C structures enumerated from derivative superstructures with volumes times the volume of the base 4 site hexagonal cell. It turns out that the experimental ground state structure of V 2 phase, does not correspond to the lowest energy structure in our DFT calculations. Instead, there are several ordered V 2 C s tructures that are ~0 2 meV/atom lower in energy, see Fig. 3 2. Note e nergies in Fig. 3 2 group Pbcn ). The low energy structures share the feature that no carbon atoms sit directly atop one another in adjacent plane s perpendicular to the c axis, as is the case phases. The structures with the lowest energies for each supercell volume possess linear chains of carbon atoms perpendicular to the c axis. Such chains are arranged parallel to one another forming a pattern in which the length of the motif increases with the volume of the structure. The simplest of these structures has a motif length of two, belongs to the orthorhombic space group Pnnm and has a volume twice that of the base hexagonal cell, see Fig. 3 3 (a). The energy of the Pnnm structure is
43 L ocal environments, out to 2 nd nearest neighbors ( NNs ) of the larger structures with longer motifs ar e very similar to that of the Pnnm structure. Therefore, we expect the self diffusion behavior of these structures to phases of V 2 C, the low energy Pnnm structure is included in our subsequent calculations an d discussion. The unit cells for these three structures are shown in Fig. 3 3. In the figure, vanadium atoms are shown as large (red) spheres and carbon atoms as small (brown) spheres. phase is unique in that the density of carbon planes p erpendicular to the c axis alternates between a low and high value, whereas this density is constant for the two orthorhombic structures. Figure 3 3. Unit cells of three ordered V 2 C structures including (a) the low energy orthorhombic structure belongin g to space group Pnnm Zero temperature lattice constants for the low energy orthorhombic structure and phases of V 2 C calculated using the GGA PBE functional are given in Table 3 1. phase are also shown. To the best of our knowledge, no phase. The equilibrium lattice parameters for phases of V 2 C calculated using the GGA PBE functional are in excellent
44 agreement with experiment, while those obtained using the LDA functional underestimate experiment by ~2%. Therefore, the GGA PBE functional is used for all calculations presented in the remainder of this chapter 3.3.2 Electronic and Elastic Properties of V 2 C 2 C at zero temperature are shown in Fig. 3 4. Band energie s are given relative to the Fermi energy which is set to zero The absence of a band gap is in agreement with the experimental fact that early transition metal carbides are good electrical conductors  Both e phase and the Pnnm structure projected DOS calculations show that states near the Fermi level are predominantly associated with vanadium atoms. Since all three ordered V 2 C structures share an almost identical vanadium sublattice, it is not surprising that the ir electronic structures are similar. Small differences can be attributed to the differing arrangements of carbon atoms. Figure 3 4. 2 C from DFT calculations.
45 Table 3 2. Zero temperature elastic constants and elastic moduli for the three ordered V 2 C structures. Elastic constants for the three different V 2 C structures are given in Table 3 2 There are no experimental measurements of single crystal elastic constants for either  and  are within 2% of those calculated in this work, while those of Ref.  are significantly higher overall. This difference may be due to the use of norm conserving pseudopotentials in Ref.  while our calculations and those of Ref.  use PAW pseudopotentials. It is not clear what type of pseudopotentials were used in the Space group (phase) Pnnm Pbcn P 3 1m ( ) 406 384 381 450 383 418 425 415 410 493 414 418 434 404 393 452 400 431 150 188 189 205 189 176 116 182 181 207 182 148 161 124 122 146 120 148 85 107 107 122 110 131 153 128 125 143 130 131 95 136 131 161 135 121 235 243 250 279 242 246 119 120 118 140 121 130 305 309 303 359 311 331 0.284 0.288 0.290 0.290 0.286 0.277 Reference This work (VASP PBE) This work (VASP PBE) VASP PBE Ref.  CASTEP PBE Ref.  CASTEP PBE Ref.  This work (VASP PBE)
46 calculations of Ref.  phase. All values a re obtained for zero temperature. Figure 3 5. E for (a) low energy orthorhombic structure belonging to space group Pnnm 2 C phase of V 2 C The bulk modulus K and shear modulus G for the aggregate V 2 C crystal a re estimated from the elastic constants using the Voigt Reuss Hill approximation  The Young s modulus and Poisson ratio a re then calculated from the well k nown relations: and Elastic moduli from our DFT calculations for all three V 2 C structures as well as from other theoretical works are given in Table 3 2. Overall the three V 2 C structures exhibit elastic moduli that differ by < 10%, which indicates their average elastic response is similar and only weakly influenced by carbon ordering. However, the two orthorhombic structures exhibit greater anisotropy of E G and v phase, which has a nearly isotropic elastic response. To illustrate this, a 3D plot of E is shown in Fig. 3 5 for each of the three V 2 C structures, where the distance from the center to the surface of each shape corresponds to E in that specific direction. More spherical shapes correspond to more isotropic elastic behavior. It is clear from the figure that E phase is nearly
47 isotropic, whereas the other two structures exhibit more anisotropy. Similar behavior is observed for G and v (plots not shown) Figure 3 6. Theoretica l phonon DOS for the (a) low energy orthorhombic structure belonging to space group Pnnm 2 phase of V 2 C. Figure 3 6 shows the theoretical phonon DOS for all three of the ordered V 2 C structures c alculated using the relaxed lattice parameters in Table 3 1. Each of the phonon DOS decays to zero as the frequency approaches zero indicating no unstable modes. Therefore, the three ordered V 2 C structures are mechanically stable. Zero point vibrational e nergy corrections to the total energy are given by where is the normalized phonon DOS  The zero point correction to the cohesive energies for each of the three V 2 C structures are given in Table 3 1. The three structures give values of E ZP that are the same to within 1 meV/atom. Therefore, the Pnnm structure remains the lowest energy V 2 C structure, though the energy difference ( 3 2 )
48 3.3.3 Vacancy F ormation E nergies As a first step towards investigating the self diffusion behavior of V 2 C, we calculated the formation energies of intrinsic defects most likely to aid in diffusion. The formation energy at zero temperature for a neutral defect X is given by [179 182] ( 3 3 ) where is the total energy of the cell with the defect, is the total energy of the perfect cell, is the number of atoms of chemical species j that are removed ( ) or added ( ) to the cell and is the energ y of the atomic reservoir to which the atom is removed or added. In this work, is taken to be the co hesive energy per atom of either bcc vanadium or carbon graphite. Figure 3 7. Calculated vacancy formation en ergies as a function of the number of atoms in the supercell for the (a) low energy structure belonging to space group Pnnm 2 phase of V 2 C. Since the self diffusion of vanadium and carbon are assumed to occur separatel y on their respective sublattices, the only relevant defects to consider are vacancies. Figure 3 7 shows the convergence of the carbon and vanadium vacancy formation
49 energies with respect to supercell size for the each of the three V 2 C structures. The vaca ncy formation energies listed in Table 3 3 a re found by extrapolating linear fits of the calculated formation energies as phase, we denote the carbon lying in the high density plane as C and the carbon lying in the low density plane as C. In the limit and differ by ~0.03 eV. Table 3 3. Vacancy formation energies and vibrational formation entropies for the three ordered V 2 C structures calc ulated using DFT Creation of a vanadium vacancy requires breaking bonds with 12 NN vanadium atoms, whereas creation of a carbon vacancy requires breaking bonds with only 6 NN vanadium atoms. Therefore, more energy is required to create a vanadium vacancy than a carbon vaca ncy and Carbon vacancy formation energies for the three different ordered structures are within 0.1 eV of each other, which is reflective of their nearly identical vanadium sublattices. Vanadium vacancy formation energies are also simi lar across structures which indicates that the carbon arrangement has only a small effect on the overall bonding strength. To the best of our knowledge, no other theoretical calculations or experimental measurements of the vacancy formation energies have been performed for V 2 C or any of the other vanadium carbides. Space group (phase) Pnnm Pbcn P 3 1m ( ) (eV) 2.42 2.48 2.43 (eV) 2.40 (eV) 4.00 4.08 4.24 ( k ) 5.00 5.53 4.97
50 3.3.4 Energy Barriers for Self diffusion 18.104.22.168 V anadium s elf d iffusion Ignoring the small distortions in the two orthorhombic structures, the vanadium atoms in V 2 C form an hcp sublattice. It is assumed that vanadium self diffusion occurs exclusively on the vanadium sublattice by a standard vacancy assisted mechanism. T his assumption is based on experimental evidence that the metal atom diffusivity in transition metal carbides i s independent of carbon concentration. Considering only the first NN shell, a vanadium may either migrate within its own basal plane (jump A) or into one of the two adjacent planes (jump B) as shown in Fig. 3 8. In the figure v anadium atoms are shown as l arge (red) spheres, while octahedral sites are shown as small (brown) spheres. In V 2 C, only half of the octahedral sites are occupied by carbon atoms The in plane and out of plane self diffusion coefficients for vacancy mediated diffusion are give n by where is the concentration of thermal vacancies on the vanadium sublattice, are the effective frequencies, and are the vacancy migration energies for jumps of type A and B, respectively, and are the in plane and out of plane squared displacements and s are the correlation factors. The subsc ripts i and j denote the index of the migrating atom. In the expression for the sum in i is over the six NN vanadium atoms within the basal plane of the vacancy and the sum in j is over the six NN ( 3 4 )
51 vanadium atoms in the two adjacent basal planes. In th e expression for the sum is taken over only the six NN vanadium atoms in the two adjacent basal planes; jumps within a basal plane do not contri b ute to Figure 3 8. In plane (A) and out of plane (B) jumps for vanadium to a nearest neighbor site within the vanadium sublattice (a) without octahedral sites shown and ( b) with octahedral sites shown. The concentration of vacancies on the vanadium sublattice follows an Arrhenius dependence Contributions from lattice vibrations to the free energy are treated in the harmonic approximation. The vibrational entropy of a harmonic solid of N atoms is given by where is the normalized phonon DOS at zero temperature  The expression for is analogous to Eq. 3 3 for and requires cal culatin g for three system s: perfect V 2 C, the same system with a single vanadium vacancy and bcc vanadium. The supercell sizes used in the force constant calculations are (96 ( 3 5 )
52 atoms) for the Pnnm structure, (72 phase A conventional supercell is used in the force constant calculation for bcc centered mesh of phonon wave vectors is used to sample the dynamical matrices, which are solved to get for each system. The integration in Eq. 3 5 is then carried out to get the vibrational entropy Calculated values of for the three structures are given in Table 3 3. The free energy of vacancy formation considering only harmonic contr ibutions from phonons is given by Migration energy is defined as where is the energy of the saddle point along the MEP between the two end states and is the energy of the initial state. For the vacancy mechanism considered here, the initial state corresponds to a supercell with a single vanadium vacancy, while the final state corresponds to one of the 12 NN vanadium atoms moved to the vacancy site. CI NEB calculations are performed on a set of in terpolated structures for each of the 12 initial and final state pairs. We find that migration energies converge rapidly with respect to supercell size. A 71 phase and 95 atom supercells for the two orthorhombic structures give converged to within 0.05 eV. MEPs for the six in plane and six out of plane migration paths for the of V 2 C are shown in Fig. 3 9( a) and 3 9( b), respectively. Three distinct groups of MEPs are observed for the out of plane B jumps shown in Fig. 3 9( b). phase and the Pnnm structure plane migration energies are found to vary by 0.2 eV, while out of plane migration energies vary by as much as
53 2.6 eV. The large variation in is attributed to the vanadium atom having to pass through an intervening plane of carbon atoms during its jump to an adjacent basal plane; see the Discussion for more details. The other two V 2 C structures show similar variation in whereas the variation in for the Pnnm structure (3.6 eV) is ~1 eV attributed to differences in carbon atom arrangements among the structures a nd will be expanded upon in the Discussion. Table 3 4 gives calculated values for and for each of the six possible jumps of types A and B. Figure 3 9. Theoret ical MEP s for vanadium vacancy migration for (a) in plane (jump A) and (b) out of plane (jump B) migration for CI NEB calculations The effective frequencies for an N atom system are calculated from where the and are the normal mode frequencies for the initial and transition states, respectively [29,30,37] The product over the transition state frequencies is reduced by one relative to the initial state due to an imaginary mode along the migration ( 3 6 )
54 direction. We approximate the normal mode frequencies for the entire system by those of the migrating atom only. This is the Einstein (independe nt oscillator) approximation. Eq. 3 6 simplifies greatly as the migrating atom contributes only three normal modes to the initial state and two modes to the transition state; the third mode of the transition state is imaginary and does not enter Eq. 3 6. T he are calculated for each of the 12 MEPs determined by the CI NEB calculations. We use a supercell with the same dimensions as was used in the calculation of and The calculated ranges for for both A and B jumps are given in Tabl e 3 4. T able 3 4. Migration energies and and effective frequencies for the vacancy assisted vanadium self diffusion mechanism for the three ordered V 2 C structures calculated using DFT Space group (phase) Pnnm Pbcn P 3 1m ) ) (eV) 2.69 2.64 2.76 2.69 2.64 2.76 2.85 2.71 2.85 2.85 2.72 2.85 2.85 2.79 2.87 2.85 2.82 2.87 (eV) 1.24 1.46 1.59 2.09 1.98 1.59 2.09 1.99 1.69 2.10 1.99 2.10 2.10 2.01 2.10 4.84 4.02 4.34 (THz) 9.3 10.0 8.3 10.5 8.2 10.9 (THz) 5.6 12.8 5.7 10.9 5.4 10.6
55 The correlation factors f for vacancy assisted self diffusion in an hcp lattice have been calculated elsewhere  and are a fu nction of the ratios of in plane and out of plane jump frequencies and respectively. Table 3 5 gives correlation factors over a temperature range of 1000 3000 K calculated using the average jump fr equencies and for vanadium self phase of V 2 C. The differences between f Bx and f Bz are negligible for the range of jump frequencies considered. The correlation factors for the other two structures are calculated in the same way and we obtain values nearly identical to those given in Table 3 5 phase. This is expected as the av erage in plane and out of plane migration barriers and effective frequencies are similar for all three V 2 C structures. Table 3 5. Correlation factors for vacancy assisted self diffusion in an hcp lattice taken from Ref.  calculated using values of obtained for vanadium self phase of V 2 C (K) ( ) 1000 0.0050 0.9214 0.6475 1200 0.0123 0.9197 0.6490 1400 0.0233 0.9171 0.6514 1600 0.0378 0.9138 0.6546 1800 0.0549 0.9097 0.6584 2000 0.0740 0.9052 0.6628 2200 0.0946 0.9004 0.6674 2400 0.1160 0.8955 0.6722 2600 0.1379 0.8906 0.6770 2800 0.1598 0.8857 0.6817 3000 0.1817 0.8810 0.6863
56 In plane and out of plane self diffusion coefficients and for vanadium are calculated by plugging in t he values of , , and listed in Tables 3 3 and 3 4 into Eq. 3 4. The self diffusion coeffic ients are then fit to an Arrhenius expression of the form Activation energies and diffusion prefactors from these fits are given in Table 3 6 for each of the three ordered V 2 C structures. Experimental values for and for meta l atom self diffusion in cubic NbC x and TiC x are also given for comparison. Overall, there is good agreement between the calculated activation energies and that measured in self diffusion experiments for Nb in NbC x while the experimental activation energy for self diffusion of Ti in TiC x is significantly higher. Calculated diffusion prefactors underestimate experimental results by an order of magnitude. This is not surprising as first principles studies of self diffusion in hcp metals [33,37] also observed a similar underestimation in their calculated diffusion prefactors. Table 3 6. Activation energies and diffusion prefactors for vacancy assisted vanadium self diffusion for the three ordered V 2 C structures calculated using DFT 22.214.171.124 C arbon s elf diffusion In each of the V 2 C structures, carbon atoms fill half of the available octahedral interstitial sites, while the other half are empty. Because the concentration of such This work ( Pnnm ) This work This work phase) Nb in NbC x (Expt. Ref.  ) Ti in TiC x (Expt. Ref.  ) (kcal/mol) 121.4 129.3 139.6 (kcal/mol) 121.2 128.9 139.5 0.089 0.223 0.238 0.383 0.812 0.852
57 structural vacancies dwarfs the concentration of thermal vacancies, we expect carbon self diffusion to occur through migration of carbon atoms to empty NN octahedral sites. Together the octahedral interstitial sites form a simple hexagonal lattice (ignoring the small distortions for the two orthorhombic structures). Coordination polyhedra for the octahe dral and tetragonal interstitial sites are shown in Fig. 3 10 C arb on atoms are surrounded by brown octahedra, while empty octahedral sites are surrounded by blue octahedra. Tetrahedral sites are surrounded by white tetrahedra. Octahedral sites share faces with neighboring octahedral sites in adjacent basal planes, but share edges with neighboring octahedral sites within the same basal plane. Tetrahedral sites share faces with neighboring octahedral sites within the same basal plane. Figure 3 10 C oordin ation polyhedra for the (a) octahedral and (b ) tetrahedral interstitial sites in V 2 C. It has been argued that in close packed structures migrating atoms will pass through faces of connecting polyhedra rather than edges due to the larger effective provided by the faces [185,186] In light of this, an octahedral tetrahedral octahedral (O T O) mechanism was proposed [20,24] in which the migrating interstitial atom temporarily occupies a NN tetrahedral site before advancing to an empty
58 octahedral site within the same basal plane. However, to the best of our knowledge, no direct evidence for the carbon self diffusion mechanism has been given for V 2 C. Therefore, we consider a direct octahedral octahedral (O O) mechanism in addition to the O T O mechanism in our DFT calculations. Between basal planes the octahedral sites share faces and so a direct O O jump mechanism is almost certain. In total, we consider two O O and one O T O jumps: an O O jump within the same basal plane (jump A), an O O jump from one basal pl ane to another (jump B) and an O T O jump in which the two octahedral sites share the same basal plane (jump C). These jumps are shown in Fig. 3 11 In the figure, v anadium atoms are shown as large (red) spheres, octahedral sites as medium sized (brown) spheres and tetrahedral sites as small (white) spheres. Figure 3 11 Octahedral octahedral, labeled A and B, and octahedral tetrahedral octahedral, labeled C, site jumps considered as likely migration paths for carbon self diffusion in V 2 C (a) without a nd (b) with the vanadium atoms shown. Based on our DFT calculations, the energy cost of a carbon atom occupying a tetrahedral site is > 1.5 eV more than the migration energy for the direct in plane O O
59 mechanism (jump A). This is true for all three of the ordered V 2 C structures considered in this work. Therefore, the contribution of the O T O mechanism (jump C) to the self diffusion coefficient of carbon will be negligible in comparison to the contribution from the direct O O mechanism (jump A). Ignoring t he contribution from O T O jumps, the in plane and out of plane self diffusion coefficients for interstitial diffusion of carbon are given by where and are the fractions of vacant octahedral sites surrounding the migrating carbon atom within the same basal plane and in the adjacent basal planes, respectively Effective frequencies are given by migration energies for jumps A and B are given by and and f the sum is taken over all possible jumps of type A. Similarly, the sum in is taken over all possible jumps of type B. Jumps of type B do not contribute to while jumps of type A do not contribute to phase, for the low density carbon plane and for the high density plane. for the two orthorhombic structures and for all three structures. Every carbon atom is surrounded by eight NN octahedral sites: six sites in the same basal plane and one site in each of the two adjacent basal planes. For all three structures, the two sites in the adjacent planes are empty, thereby contributing two jumps of type B. For the two orthorhombic structures, four of the six in plane sites are empty, while two are occupied by carbon atoms. Thus, there are four possible A jumps ( 3 7 )
60 phase, the situation is complicated by t he fact that there are low and high density carbon planes. In the low density plane, all six in plane sites are empty, whereas, in the high density plane, only three in plane sites phase, there are six possible A jumps for the low density carbon plane and three possible A jumps for the high density plane. The converged phase are shown in Fig. 3 12 The MEPs for the orthorhombic structures (not shown) are similar albeit with only a single uniqu e carbon plane. As in the case of the CI NEB calculations for vanadium migration, the saddle points for carbon migration converge rapidly with respect to supercell size. Migration energies for a 72 phase and 96 atom supercells for the two orthorhombic structures are converged to within 0.05 eV. Figure 3 12 Theoretical MEPs for carbon interstitial migration (a) in plane (jump A) and (b) out of phase from CI NEB calculations The converged migration ene rgies and for each possible jump of that type and for each structure are given in Table 3 7. There is a variation in of 0.2 eV are equal for the Pnnm structure and in like phase. These differences can be attributed to the distribution
61 of carbon atoms among the six in plane NN octahedral sites surrounding the migrating carbon atom. In the Pnnm phase, the distribution o f carbon atoms is and for carbon self diffusion are also given in Table 3 7 for each jump type. The for carbon are calculated in the same manner and u sing the same supercell sizes used in the calculation of for vanadium self diffusion. Table 3 7. Migration energies and and effective frequencies for carbon self diffusion for each of the three ordered V 2 C structures calculated usi ng DFT The vanadium and carbon atoms do not share a common sublattice: one is hcp while the other is simple hexagonal. This contrasts with fcc interstitial compounds, where the octahedral interstitial sites also form an fcc lattice. In the fcc case, the Space group (phase) Pnnm Pbcn P 3 1m ( ) (eV) 4.09 4.13 4.44 (low density) 4.10 (high density) 4.09 4.18 4.44 (low density) 4.10 (high density) 4.09 4.18 4.44 (low density) 4.10 (high density) 4.09 4.37 4.44 (low density) 4.44 (low density) 4.44 (low density) (eV) 3.57 3.32 3.30 (low density) 3.62 (high density) (THz) 19.8 19.9 20.0 23.4 22.6 22.8 (low density) 21.7 21.8 (high density) (THz) 19.7 19.8 21.6 21.7 20.5 20.6 (low density) 22.0 (high density)
62 correlation factor for interstitial self diffusion will lie between that for pure vacancy assisted diffusion and unity. Because there is not a simple one to one correspondence for hcp structures we instead set the correlatio n factors for carbon self diffusion to unity. The self diffusion coefficients given in Eq. 3 7 are calculated using the results in Table 3 7 and then fit to The activation energies and diffusion prefactors from these fits are given in Table 3 8. Overall, the calculated activation energies for all three ordered V 2 C structures agree well with the experimental results for self diffusion of C in disordered cubic V 6 C 5 and C in cubic NbC 0.868 However, our calculated diffusion prefactors underestimate experiment by two orders of magnitude. Table 3 8. Activation energies and diffusion prefactors for carbon self diffusion for each of the three ordered V 2 C structures calculated using DFT 3.4 Discussion W hen a vanadium atom migrates from one basal plane to an adjacent one, it must pass through an intervening plane of octahedral sites. The black arrow labeled B in Fig. 3 8 (b) shows the direction of the out of plane jump. At the transition state, the vanadiu m atom lies between two octahedral sites. This leads to one of three scenarios. This work ( Pnnm ) This work This work phase) C in V 6 C 5 (Expt. Ref.  ) C in NbC 0.868 (Expt. Ref.  ) (kJ/mol) 91.5 96.3 102.4 (low density) 94.5 (high density) (kJ/mol) 82.9 76.6 76.1 (low density) 83.4 (high density) 0.016 0.015 0.028 (low density) 0.014 (high density) 0.009 0.010 0.009 (low density) 0.010 (high density)
63 Either both octahedral sites are occupied by carbon atoms, both sites are vacant, or only one site is occupied. Due to the proximity of the octahedral sites to the migrating v anadium atom at the transition state, the height of the migration barrier depends strongly upon the distribution of carbon atoms at these sites. In the two orthorhombic structures, where the intervening plane of carbon atoms always has the same density, th e highest migration barrier occurs when both sites are occupied, while the lowest barrier occurs when both sites are vacant. The third scenario, in which only one site is occupied, leads to an intermediate barrier. Evidence for this behavior is given by th e three groups of MEPs shown in Fig. 3 9( Pnnm structure exhibits the same phase is complicated because of the alternating high and low density carbon planes. However, the highest out of plane vanadium phase still corresponds to both octahedral sites being occupied by carbon atoms. Despite the two o rthorhombic structures exhibiting the same grouping of MEPs for out of plane vanadium jumps, we found that that the barrier heights for jumps within a group are different for each structure. For example, the highest and lowest out of plane vanadium migrati respectively, while, for the Pnnm structure, they are 4.84 eV and 1.24 eV. To understand this behavior, one must also consider the arrangement of other nearby carbon atoms in addition to those closest to the vanadium migration path. Carbon atoms in the Pnnm structure form linear chains that lie parallel to basal planes of vanadium
64 chains is highlighted in Fig. 3 13 In the figure, l arge (red) spheres are vanadium atoms, small dark (brown) spheres are carbon atoms and small light (white) spheres are empty octahedral sites. Figure 3 13 V 2 C and (b) the low energy orthorhombic structure described by space group Pnnm A vanadium atom migrating into an adjacent basal plane in which the two intervening octahedral sites are occupied will displace both carbon atoms away from their ideal sites. O n one hand, for the Pnnm structure, the displacements are along the direction of the carbon chain. This leads to a great deal of resistance from the strong carbon carbon bonds and the corresponding migration barrier for the vanadium atom is high. On the ot closest to the migrating vanadium atom are displaced towards empty octahedral sites. This leads to less resistance from other carbon atoms in the chain and correspondingly a smaller migra Pnnm structure. The lowest out of plane migration barriers occur when both octahedral sites closest to the migrating vanadium atom are empty. At the transition state, the closest carbon atoms to the migrating vanad ium atom are displaced towards the partially
65 vacated vanadium site. In the Pnnm structure, the displacements are nearly perpendicular to the direction of the carbon chains, so there is minimal resistance from other carbon atoms within the chain. By contras these nearby carbon atoms tend to be directed more towards other neighboring carbon atoms within the zigzag chains. This leads to more resistance due to the strong carbon carbon bonds and the corresponding migration the Pnnm structure. Because the lowest migration barrier dominates the self diffusion coefficient, the calculated activation energy for vanadium self diffusion is ~8 kJ/mol less for the Pnnm structure than for th Of the three V 2 phase has the highest activation energy for vanadium self Pnnm structure, see Table 3 6. phase is largely due to a ~0.2 eV difference in calculated values of see Table 3 3 The calculated self diffusion coefficients for vacancy assisted self diffusion of vanadium for the three ordered V 2 C structures are plotted in Fig. 3 14 Experimental self diffusion data for Nb in NbC x  and Ti in TiC x  are also included in the plot for comparison. Overall, the calculated self diffusion coefficients are in better agreement with the self diffusion data for Nb in NbC x than for Ti in TiC x which is expected, as Nb and V belon g to the same group in the periodic table.
66 Figure 3 14 Calculated vanadium self diffusion coefficients for the three ordered V 2 C structures and experimental self diffusion coefficient s for Nb in NbC x  and Ti in TiC x  The calculated self diffusion coefficients for interstitial diff usion of carbon in the three ordered V 2 C structures are shown in Fig. 3 15 along with experimental self diffusion data for C in cubic NbC 0.868  and C in cubic V 6 C 5  The self diffusion coefficients and phase are averaged over the low and high density carbon planes. Note that the experimental data was obtained for self diffusion in disordered structures, not the ordered structures considered in this work. Overall the calculated activation energies in Table 3 8 agree well with the experimental result for C in V 6 C 5 However, the calculated diffusion prefactors underestimate the experimental values by two orders of magnitude, which results in bad overall agreement with the self diffusion coefficients, which is clear from Fig. 3 15 It is unclear why the underestimation is so sever e.
67 Figure 3 15 Calculated carbon self diffusion coefficients for the three ordered V 2 C structures and experimental self diffusion coefficients for C in NbC 0.868  and C in V 6 C 5  Comparison of the vanadium and carbon self diffusion coefficients shows that carbon self diffusion is se veral orders of magnitude higher than vanadium self diffusion. This is in agreement with self diffusion experiments on several cubic transition metal carbides in which the carbon diffusivity was also found to be several orders of higher than the metal atom diffusivity [16,20,24,26 28,158] Because we assume carbon diffuses via an interstitial mechanism, the carbon vacancy formation energy does not enter the expression for the self diffusion coefficient. In contrast, vanadium self diffusion is mediated by thermal vacancies and the added energy cost to produce these vacancies results in a much lower self diffusion coefficient relative to carbon.
68 3.5 Summary By enumerating all possible ordered V 2 C structures up to 8 times the volume of the base hcp unit cell, our calculations reveal several ordered structures with lower en phases. The simplest of these low orthorhombic space group Pnnm The three ordered structures all exhibit similar electronic and elastic properties due to their nearly identical vanadium sublattices with variations that can be attributed to their different carbon arrangements. Of the three ordered V 2 phase has the highest ac tivation energy for vanadium self Pnnm structure. Differences in activation energies among the three ordered structures are the result of a combination of differences in vanadium vacancy formation energies an d migration barrier heights. Calculated vanadium self diffusion coefficients agree well with experimental results for Nb self diffusion in single crystal samples of disordered cubic NbC x  The calculated carbon self diffusion coefficients are several orders of magnitude larger than those calculated for the metal atom, which is in agreement with diffusion experiments for the cubic transition metal carbides [16,20,24,26 28,158] Theoretical activation energies for carbon self diffusion in ordered V 2 C agree well with experimental results for carbon self diffusion in disordered cubic V 6 C 5 [ 24] while diffusion prefactors are severely underestimated.
69 CHAPTER 4 ENERGETICS AND KINETICS OF METAL IMPURITIES IN VANADIUM CARBIDE 4.1 Background Diffusion experiments on transition metal carbides are difficult due to a number of factors, including the difficul ty of producing the high purity single crystal samples needed for accurate diffusion measurements  While self diffusion behavior of carbon interstitials has been investigated experimentally in NbC [20,23] ZrC  TiC [16,17,22] V 6 C 5  WC  Mo 2 C  and Fe 3 C  self diffusion measurements for the metal atom are quite rare and results are known only for NbC x  and TiC x  Furthermore, experiments have yet to address the diffu sion of impurities in any of the transition metal carbides, including V 2 C. From a theoretical perspective, the energetics of intrinsic defects and self diffusion behavior in ordered phases of V 2 C have been addressed in Chapter 3 However, no work has been done to address the energetics of impurities in V 2 C or their diffusion behavior. In this chapter we seek to investigate the energetics of metallic impurities in the dilute limit and their migration within V 2 C using DFT calculations. As was discussed in Ch apter 3, this study is motivated by the potential use of V 2 C as a diffusion coating to help mitigate fuel cladding chemical interaction in nuclear fuel rods. Calculations are performed for the ordered orthorhombic phase V 2 C as (1) application temperature s for the diffusion barrier coating are expected to be below the order disorder transition temperature (2) it is the confirmed ground state phase for V 2 C and (3) results of Chapter 3 indicate that carbon ordering does not have a significant effect on the self diffusion behavior Impurity species considered in this work include Fe and Ni as major
70 components of the cladding, Ce and Nd as fission products and U as a major component of the nuclear fuel. As a first step, the formation energies for impurities at both the vanadium substitutional site and octahedral interstitial site are calculated to determine the site preference for each chemical species. All im purities are found to prefer the vanadium substitutional site. Assuming a standard vacancy assisted mechanism for diffusion of impurities we then calculate the impurity vacancy binding energies. Large impurities, U, Nd and Ce are found to bind much more strongly to neighboring vanadium vacancies than the smaller impurities, Fe and Ni. Relaxed impurity vacancy structures show the large impurities form clusters with neighboring vacancies, so a standard vacancy assisted diffusion mechanism is deemed unlikely M igration barriers are calculated for the two small impurities, Fe and Ni, as they do not aggregate with neighboring vacancies. Comparison with migration barriers for V from Chapter 3 show that migration is easier for Ni, whereas it is harder for Fe. 4.2 Simulation Methods Calculations a re performed using the plane wave DFT code VASP [131,159] PAW potentials a re used to treat the core electrons [160,161] Electron exchange correlation i s treated by the PBE  formulation of the GGA which has been shown in Chapter 3 to outperform the LDA in calculations on V 2 C. V alence electron configurations for the reference atoms a re the same as in Chapter 3 with the addition of 3 d 6 4 s 2 for Fe, 3 d 8 4 s 2 for Ni, 5 s 2 5 p 6 5 d 1 6 s 2 for Ce, 5 s 2 5 p 6 5 d 1 6 s 2 for Nd and 6 s 2 6 p 6 5 f 2 7 s 2 6 d 2 for U DFT is known to have difficulties in treating partially filled f orbitals. Therefore, in our calculations t he 4 f electrons for Ce and Nd a re treated as part of the frozen core in the pseudopotentials as was done in previous DFT calculation of
71 rare earth im purities in Mg  There is no available PAW potential for U that includes the 5 f electrons in the frozen core. Therefore, we apply a Hubbard U correction to the 5 f orbitals in U using the Dudarev formulation  A n effec tive Hubbard parameter U eff = U J = 1.24 eV i s used, which was shown to be the optimal value for calculations involving U metal and U Zr alloys  Tolerances on the electronic self consistency cycle, interatomic forces and cel l stress tensor are the same as in Chapter 3. Partial orbital occupancies a re treated using the method of Methfessel Paxton  (2 nd order) w ith a smearing width of 0.1 eV. Table 4 1 Theoretical equilibrium lattice constants, cohesive energies and magnetic moments for ground state phases of V, Fe, Ni, Ce, Nd and U Deltas signify percent differences from experiment  Note: ortho orthorhombic; dhcp double hexagonal close packed Reference chemical potentials used in the calculation of defect formation energies a re calculated from the total energies of the relaxed ground state crystal structures for each element. Spin polarized calculations a re performed for the ferromagnetic elements Fe and Ni. Ground state structures, equilibrium lattice parameters, cohesive ener gies and magnetic moments (whe re necessary ) are given in V Fe Ni Ce Nd U Structure bcc bcc fcc fcc dhcp ortho a () 2.997 2.832 3.518 4.737 3.705 2.803 0.9% 1.2% 0.2% 2.3% +1.3% 1.8% b () 5.837 +0.6% c () 11.926 4.911 +1.1% 0.9% E 0 (eV/atom) 8.992 5.468 8.238 5.926 4.712 11.136 B ) 2.20 0.65 +0.5% +2.3%
72 Table 4 1 along with percent age differences from experiment  A plane wave energy cutoff of 600 eV and 6000 k points per reciprocal atom a re found to give total energies converged to within a few meV. Overall the agreement between the calculated and experimental lattice constants is good with a maximum difference of a few percent. Magnetic moments fo r bcc Fe and fcc Ni are also within a few percent of experiment. Figure 4 1. Unit cell (solid lines) of low temperature ordered orthorhombic phase V 2 C. The vanadium sites of the V 2 C phase constitute a nearly hcp sublattice with carbon atoms occupying exactly half of the availabl e octahedral interstitial sites. The unit cell of V 2 C is shown in Fig. 4 1 with additional vanadium atoms included to emphasize the hcp like vanadium sublattice. In the figure, vanadium atoms are shown as large (red) spheres and carbon atoms are shown as small (brown) spheres. The octahedral site sublattice on which the carbon atoms reside is simple hexagonal with basal planes layered between the vanadium basal planes. Note that the close packed direction is along which means the basal planes lie in (100) lattice planes. Table 4 2 gives the cell dimensions, number of atoms, lattice parameters and k point grids for the
73 supercells used in the defect calculations. Supercells up to 324 atoms are used to test for conver gence of defect formation energies. The minimum distance between defects is given by the smallest cell dimension, e.g. a single defect in a 222 supercell has its nearest periodic image 9.10 away. Table 4 2 Supercell dimensions and calculation paramet ers used in defect calculations for V 2 C MEPs for impurity migration are determined using the CI NEB method [192,193] The CI NEB is preferred over the standard NEB method as one of the images along the MEP is guaranteed to converge to the saddle point, thus allowing easy and accurate calculation of migration barriers. The spring constant between images is set to 5 eV/ 2 Convergence of the MEP is tested with respect to the number of interpolated images between the initial and final structu res for a representative migration path. It is found that the MEP, including the location of the transition state, does not change appreciably beyond five images; therefore, five images are used in all CI NEB calculations presented in this chapter Supercell size No. of atoms a () b () c () k point grid 111 12 4.55 5.73 5.03 999 122 48 4.55 11.46 10.06 955 222 96 9.10 11.46 10.06 555 322 144 13.65 11.46 10.06 355 323 216 13.65 11.46 15.10 353 333 324 13.65 17.20 15.10 333
74 4.3 Simulat ion Results 4.3.1 Impurity Site Preference The substitutional vanadium and octahedral interstitial sites are chosen as the most likely sites for an impurity to occupy. All early transition metal carbides are good electrical and thermal conductors [146,147] A metallic character is then expected for V 2 C and is confirmed by electronic band structure and DOS calculations performed in Chapter 3 (Fig. 3 4) and in Refs. [34 36 ] which show an absence of a band gap at the Fermi level. Therefore, the metallic impurities will be charge neutral. The formation energy for a neutral impurity X (= Fe, Ni, Ce, Nd and U ) on a substitutional vanadium site is given by where is the total energy of the supercell with the impurity, is the energy of the perfect supercell and and are the chemical potentials of vanadium and the impurity, respectively. The chem ical potential for each element is taken as the cohesive energy per atom of the relaxed ground state phase, see Table 4 1 All vanadium sites in V 2 C are equivalent by symmetry, so there is a single unique for each impurity species. Half of the octahedral sites are occupied by carbon atoms, while the remaining half are free to accept an impurity. The formation energy for an impurity at an interstitial site is given by where only the chemical potential of the impurity must be accounted for as the interstitial site is not previously occupied. Again, is the cohesive energy per atom of ( 4 1) ( 4 2)
75 the ground state phase of impurity X As is the case for the vanadium sites, all empty interstitial sites are equivalent by s ymmetry, so there is a single unique for each impurity species. Using the relaxed supercells from the substitutional and interstitial defect calculations, we also calculate the impurity volume where and are the volumes of the defect and perfect supercell, respectively. gives a measure of the size of the impurity. In the case of a substitutional defect, a negative implies that the defect is smaller than original vanadium atom Figure 4 2. Convergence of (a) substitutional and (b) interstitial defect formation energies with respect to the size of the defect supercell. In each case, the zero of the relative formation energy corresponds to the size converged formation energy. Cal culations of both substitutional and interstitial impurities are performed for the supercells listed in Table 4 2 to test for convergence of the defect formation energies. Relative defect formation energies are plotted in Fig. 4 2 as a function of the numb er of atoms in the supercell, where the reference point is taken as the corresponding defect ( 4 3)
76 formation energy for the largest 333 supercell. Formation energies for both substitutional and interstitial impurities have converged to within 0.1 eV when the s upercell dimensions are 222 and to within 0.05 eV when the supercell dimensions are 322. In addition, we observe less variation in the defect formation energies for the two smaller impurities, Fe and Ni, then for the larger ones, Ce, Nd and U Table 4 3 Substitutional and interstitial defect formation energies and impurity volumes. Table 4 3 gives the substitutional and interstitial defect formation energies and impurity volumes for Fe, Ni, Ce, Nd and U calculated for the 333 supercell. No magnetic ground states were found in spin polarized calculations. For all impurity species, the substitutional site is ~2 3 times lower in energy than for the interstitial site. Therefore, we conclude that substituti on on a vanadium site is preferred for all impurity species. For both sites, the defect formation energy increases with the impurity volume. The ordering of impurity volumes is Fe < Ni < U < Nd < Ce Defect formation energies follow the same trend With out including the Hubbard correction, increases by 0.1 eV and decreases by 0.2 eV. Moreover, decreases by 0.5 and 0.8 3 for the substitutional vanadium site and interstitial site, respectively. Impurity Formation energy (eV) Impurity volume ( 3 ) substitution interstitial substitution interstitial Fe 0.77 2.83 3.7 11.6 Ni 1.38 2.89 2.4 11.9 U 2.81 6.26 10. 4 19.5 Nd 4.36 8.1 2 15.3 24.3 Ce 5.07 8.82 15.6 25.8
77 Impurity volumes for Ce, Nd and U at the substitutional site are positive indicating that they are larger than the original V atom, whereas for Fe and Ni, indicating a contraction of the cell around the impurity. Both Fe and Ni undergo minimal displacements < 0.05 from the ideal vanadium site, while the larger impurities undergo much larger displacements > 0.2 directed out of the basal plane and towards empty octahedral sites. As expected, all impurity volumes for the interstitial site are positive due to the site being p reviously unoccupied. 4.3.2 Impurity vacancy Binding Experiments have found the metal atom self diffusion in NbC x  and TiC x  to be independent of the carbon composition. Therefore, as the vanadium substitutional site is energetically favored, we expect diffusion of the impurity to occur via di rect exchange with a nearest neighbor (NN) vanadium vacancy. A substitutional impurity is surrounded by twelve NN vanadium atoms: six within its own basal plane and three in each of the two adjacent basal planes, see Fig. 4 3 In an hcp crystal, like Mg, a ll in plane impurity vacancy pairs are equivalent and all out of plane impurity vacancy pairs are equivalent [36,187,195] However, in V 2 C, the impurity does not remain at the ideal vanadium site. Moreover, the proximity to carbon interstitials is different for different impurity vacancy pairs. This leads to asymmetry among the in plane and out of plane pairs; therefore, all twelve pairs must be examined individually. The binding energy for an impurity vacancy pair is given by where is the total energy of a supercell containing a substitutional impurity and a NN vanadium vacancy and is the total energy of the same supercell with a ( 4 4)
78 single vanadium vacancy. and are the same as in Eq. 4 1. The negative sign in front of the binding energy is added by convention so that indicates favorable binding. The formation energy of an impurity vacancy pair is related to the binding energy by where is the formation energy of a single vanadium vacancy and has already been calculated for V 2 C in Chapter 3 As is a fixed value, a greater value for translates into a lower formation energy for the impurity vacancy pair. Figure 4 3. Impurity vacancy pairs in the V 2 C structure: six in plane pairs have indices 1 6 and out of plane pairs have indices 7 12. Calculated binding energies for each of the twelve impurity vacancy pairs and five impurity species are given in Table 4 4 Values are calculated for a 222 supercell. Indices correspond to the labels in Fig. 4 3, which shows the 1st NN shell around a sub stitutional vanadium impurity in V 2 C. For example, the row indexed 1 in Table 4 4 corresponds to supposing the atom labeled 1 in Fig. 4 3 was removed. Almost all binding energies are positive, indicating a favorable impurity vacancy ( 4 5)
79 interaction, sa ve for one impurity vacancy pair for Fe, which has a small negative value of 0.01 eV. T he negative binding energy may be due to the increased proximity of this impurity vacancy pair to interstitial carbon atoms relative to other pairs within the same basa l plane. T able 4 4 Impurity vacancy binding energies for each of the six in plane and six out of plane nearest neighbor sites on the vanadium sublattice. indicates favorable binding. Independent of the impurity species, is highest for the 9th out of plane impurity vacancy pair. Examination of Fig. 4 3 shows that this pair has no intervening carbon interstitials between the adjacent vanadium basal planes. For the three large impurities, Ce, Nd and U the minimum value of occurs for the 12th out of plane impurity vacancy pair, which has two interv ening carbon interstitials between the adjacent vanadium basal planes. The difference between the maximum and minimum values of is > 2 eV for Ce, Nd and U which indicates the proximity of carbon interstitials to the impurity vacancy pair can have a significant effect on binding strength Index Fe Ni U Nd Ce (eV) (in plane) 1 0.07 0.21 0.45 0.40 1.21 2 0.05 0.17 0.54 0.55 0.65 3 0.09 0.18 0.61 1.15 1.46 4 0.11 0.21 0.91 1.05 1.21 5 0.01 0.14 1.01 1.13 1.31 6 0.12 0.17 0.57 0.64 0.76 (eV) (out of plane) 7 0.15 0.26 1.31 1.43 1.65 8 0.13 0.28 1.13 1.56 1.84 9 0.22 0.60 2.84 3.53 3.77 10 0.14 0.22 1.31 1.43 1.65 11 0.12 0.18 0.85 1.32 1.60 12 0.21 0.34 0.45 0.40 0.44
80 in the case of large impurities. In contrast, for smaller impurities, Fe and Ni, the difference is much less (< 0.5 eV). For all impurity vacancy pairs, binding is much weaker for the small impurities, Fe and Ni, than for the large ones, Ce, Nd and U This can be seen clearly in Fig. 4 4(a), which plots the average in plane and out of plane impurity vacancy binding energies as a function of the impurity volume. The average distance of each impurity from the ideal subst itutional site due to relaxation is plotted in Fig. 4 4(b). The much larger binding energies in the case of Ce, Nd and U are the result of significant displacements (> 0.5 ) of the impurity towards the vacancy. In contrast, Fe and Ni remain much closer (< 0.1 ) to the ideal substitutional site. As was shown in previous calculations of impurities in Al  the lattice strain that surrounds a large impurity can be reduced through the introduction of a neighboring vacancy. Clear ly, the same mechanism is at play here. Figure 4 4. Average in plane and out of plane impurity vacancy (a) binding energies and (b) distance from the ideal vanadium site for substitutional impurities X = Fe, Ni, Ce, Nd and U in V 2 C.
81 As no experimental values exist for impurity vacancy binding energies in V 2 C, we compare our calculated energies for V 2 C with those calculated in previous works for hcp Mg [187,195] Calculated for an in plane impurity vacancy pair in V 2 C ( hcp Mg) are 0.07(0.09), 0.18(0.16), 0. 82 ( 0.02 ) and 1.10 (0.0 8 ) for Fe, Ni, Nd and Ce respectively. Note that the values given for V 2 C are averages over all six in plane pairs, as in Fig. 4 4(a). Binding energies for Fe and Ni in Mg are within a few hundredths of an eV of the average binding energies for V 2 C, whereas binding energies for the larger impur ities, Nd and Ce are much higher in V 2 C than in Mg. Although impurity vacancy binding energies were not calculated in Refs. [187,195] for U in Mg, we expect the result to be the same. We stress that the only similarity between V 2 C and Mg is that the metal atoms in both structures form an hcp lattice. The close agreement observed i n impurity vacancy binding energies for Fe and Ni is coincidental. We suspect that the large difference in binding energies between the small impurities (Fe an d Ni) and the large impurities ( Ce, Nd and U ) can be attributed to the openness of the V 2 C structure, which is the result of a partially filled interstitial sublattice. This openness in V 2 C allows large impurities room to relax, so that they can get closer and bind more strongly to a neighboring vanadium vacancy. This is confirmed by the dat a in Fig. 4 4(b), which shows that Ce, Nd and U undergo large displacements away from the ideal vanadium substitutional site. In the close packed crystal structure of Mg, there is much less room into which atoms can relax, so the binding energies for the s mall and large impurities are similar. 4.3.3 Migration Barriers for Impurity Diffusion When a vanadium vacancy is introduced into the 1st NN shell surrounding one of the large impurities, Ce, Nd and U the impurity undergoes significant relaxation towards
82 the vacancy, which results in the formation of a defect cluster. In some cases, the relaxation is so large that the impurity sits within the interstitial plane. Therefore, we expect that diffusion of Ce, Nd and U should occur via some alternate mechanism, possibly involving multiple vanadium or carbon vacancies. However, the small impurities, Fe and Ni, remain closely bound to the ideal substitutional site and, so a simple vacancy assisted mechanism may still be applicable. Figure 4 5. MEPs for migration of Fe and Ni impurities for each of the six in plane and six out of plane vanadium sites. MEPs for migration of Fe and Ni impurities to NN vanadium vacancies are calculated using the CI NEB method. All six in plane and out of plane vacancy sites must be treated as was done in the calculation of the impurity vacancy binding energies. Calculated MEPs for each of the six possible in plane and six possible out of plane jumps for Fe and Ni impurities are shown in Fig. 4 5. Energies are given relative to the in itial state and the maximum point along the MEP corresponds to the transition state. MEPs calculated for in plane jumps are more clustered than those calculated for out of plane jumps. This was also seen in Chapter 3 for the case of vanadium migration
83 in V 2 C, where it was determined that the large variation in out of plane jumps can be attributed to the intervening plane of carbon interstitials. The average in plane and out of plane migration barriers in V 2 C are 2.85 and 2.86 eV for Fe and 1.96 and 2.10 e V for Ni respectively For comparison, the average in plane and out of plane migration barriers for V in V 2 C are 2.71 and 2.24 eV respectively Therefore, migration is easier for Ni than for V, while migration is more difficul t for Fe than either Ni or V. 4.4 Summary In this chapter we performed DFT calculations to investigate the energetics and kinetics of dilute metal impurities in V 2 C. The impurities Fe, Ni, Ce, Nd and U were chosen to represent the different components of the combined fuel clad system in which V 2 C may be used as a diffusion barrier. We found that all impurity types prefer the vanadium substitutional site over the octahedral interstitial site and that the energy to incorporate the impurity increases with the impurity volume. Impurity vac ancy binding energies were calculated as a first step in determining the migration kinetics. Large impurities, U, Ce and Nd, bind more strongly to neighboring vanadium vacancies than the smaller impurities, Fe and Ni. The strong binding for the large impur ities is attributed to the openness of the V 2 C structure, which allows the lattice strain created upon insertion of the impurity to be partially relieved upon introduction of a neighboring vacancy. Because the large impurities form clusters with neighborin g vanadium vacancies, a standard vacancy assisted migration mechanism is deemed unlikely. However, the small impurities, Fe and Ni, remain close to the ideal vanadium site, so migration barriers can be calculated. Comparison with migration barriers for van adium self diffusion show migration of Ni is easier, while migration of Fe is more difficult.
84 CHAPTER 5 SHOCK COMPRESSION OF COPPER ZIRCONIUM METALLIC GLASS ES 2 5.1 Background Metallic glasses ( MGs ) have attracted a great deal of attention due to their excellent mechan ical and physical properties [2,50,197 199] Because of their good glass forming ability and low cost, Cu and Zr based MGs are of interest for engineering applications. Recently, high strain rate impact experiments have been carried out on Cu and Zr based multicomponent MGs to investigate shock compression [57,58,60,65,200] as well as spallation behavior [201,202] Hugoniot compre ssion curves of shock velocity ( U s ) versus impact /particle velocity ( U p ) exhibit bilinear behavior with kinks appearing between 10 and 20 GPa that most likely indicate a phase transition [57,58,60,65,200] Because of the short time and length sca les inherent to shock compression experiments, no one has yet directly observed the microstructural changes occurring within the MGs during shock compression, and explanations of the kinks observed in the U p U s Hugoniot curves continue to be a matter of de bate [60,200] MD simulations can provide a full atomic level picture of the shock response of the sample, while overlapping with the range of strain rates attainable in laser based shock compression experiments. Such simulations should be useful in shedding light on the high strain rate deformation behavior of MGs. A large number of MD simulation works have focused on the shock response of single crystal [82 85,88,89,92,203 205] and polycrystalline metals [81,206 208] By 2 The work described in this chapter has been published in P. Wen, B.J. Demaske, D.E. Spearot and S.R. Phillpo t, J. Mater. Sci. doi: 10.1007/s10853 017 1666 5 The first two authors contributed equally to this work.
85 contrast, only a few MD simulation works have investigated the shock response of MGs [91,93,209 211] One such study by Arman et al.  found that atoms with different local structures exhibit different resistances to shear deformation, which leads to shear localization. Chen et al.  carried out MD simulations to investigate the time dependent shock induced strength of Cu 46 Zr 46 Al 8 MG. They proposed a hard viscous composite material model to describe non equilibrium strength relaxation at the shock wave front. Jian et al.  investigated short and medium range orders, and atomic packing efficiency in Cu 46 Zr 54 MG upon shock loading. They found that the evolutions of short and medium range orders are the effects of compression combined with shock induced melting. H uang et al.  studied the ductile to brittle transition during spallation of Cu 50 Zr 50 MG. They found that the transition is controlled by the interaction between void nucleation and growth, which can be regarded as a competition between tension transformation zones and shear transformation zones ( STZs ) STZs are the fundamental unit of pla stic deformation in MGs and will be discussed in more detail in Chapter 6. Lastly, Marinier et al. used a two temperature model coupled to MD t o simulate laser ablation of Cu x Zr 100 x MGs and observed homogeneous nucleation of STZs out of the pristine glass y material  Recent experimental work and MD simulations on the binary Cu Zr glass forming system have reve aled a strong dependence of glass forming ability on the composition dependence [212 214] as well as structural and mechanical properties [215 222] Over a wide composition range ( x = 35 70), Cu x Zr 100 x MGs exhibit monotonic changes in high e nergy X ray diffraction patterns, atomic pair correlation functions, mass density and thermal stability behavior  Quasistatic compression
86 experiments on a series of Cu x Zr 100 x MGs have shown that as Cu content increa ses from Cu 50 Zr 50 to Cu 64 Zr 36 the yield strength of the MG increases while the p lasticity gradually diminishes [215,216] In this chapter the shock response of Cu x Zr 100 x ( x = 30, 50 and 70) MGs is investigated via large scale MD simulations. The response of each MG is characterized over a range of piston velocities U p = 0.125 2.5 km/s, which corresponds to shock pressures from 3 to 130 GPa. At moderate shock intensities, all three MG compositions are found to undergo an elastic plastic transition characterized by the formation of a slug gishly rising plastic wave preceded by a sharply rising and narrow elastic precursor. In contrast to crystalline metals, the yield strength of MGs under shock loading is found to depend on non shear components of the pressure tensor. Simulated U p U s Hugoni ot curves exhibit elastic and plastic branches that converge at moderate shock intensities. Hugoniot states are found to depend strongly on Cu composition with Cu 70 Zr 30 exhibiting greater resistance to plastic deformation than either Cu 50 Zr 50 or Cu 30 Zr 70 5.2 Simulation Methods 5.2.1 Interatomic Potential Verification Previous MD studies of impact or laser based shock compression or Cu Zr MGs used the Finnis Sinclair type EAM interatomic potential developed by Mendelev et al.  whereas shock compression simulations of Cu Zr Al MGs used the ternary EAM potential developed by Cheng et al.  While both the Mendelev and Cheng potentials were designed to model MG structures, neither was designed specifical ly for shock compression simulations. Therefore, it is not clear which is better suited for this purpose. Zero temperature stress strain (cold pressure) curves are a computationally inexpensive way to help determine the applicability of a given potential t o shock
87 compression simulations. A well behaved potential should exhibit a monotonic increase 2 P 2 = 0 where P is the hydrostatic pressure and is the mass density Meeting this criterion is critical for shock compression simulations as failure to do so can result in unusual behavior such as decreasing sound speed ( c 2 = P ) with compression. Figure 5 1. Zero temperature stress strain curves and sound speed curves calculated using the (a ) Cheng potential for Cu Zr ( Al )  and the (b) Mendelev p otential for Cu Zr  Cold pressure curves for uniform compression are calculated using LAMMPS  for the following crystalline phases: Cu ( fcc ), Cu 8 Zr 3 ( Pnma ), Cu 7 Zr 10 ( C 2 ca ), CuZr (B2), CuZr 2 ( I 4/ mmm ) and Zr ( hcp ). These phases are chosen to represent th e different
88 local environments atoms are likely to encounter for the different compositions of MG s Figure 5 1 is a plot of the cold pressure and sound speed curves for the Mendelev and Cheng EAM potentials. Both potentials exhibit a monotonic increase in pressure for all crystalline phases over the range 0.75 < V / V 0 < 1, where V/V 0 is the compression ratio. For the Cheng potential, all four binary phases exhibit increasing sound speed with compression (no inflection points) down to V / V 0 = 0.78, while the t wo unary phases, fcc Cu and hcp Zr, exhibit decreasing sound speed with compression over the brief intervals 0.88 < V / V 0 < 0.90 and 0.84 < V / V 0 < 0.85, respectively. For the Mendelev potential, the cold pressure curves are generally not as well behaved. All four binary phases and fcc Cu exhibit decreasing sound speed with compression over a large fraction of the range 0.75 < V / V 0 < 1. The only ex ception is hcp Zr, which has similar behavior for both potentials Compression waves with amplitudes within the range where sound speed decreases with compression will not form shocks, but will instead spread out during propagation. Moreover, rarefaction w aves in this range will form shocks. Such behavior is atypical [54,225] If instead sound speed increases with compression, as is the case in m ost materials, then an initially smooth compression wave will steepen as it propagates until it forms a shock. The large regions of decreasing sound speed with compression exhibited by the Mendelev potential for both unary and binary crystalline phases of the Cu Zr system make it unattractive for shock simulation. In contrast, the Cheng potential exhibits decreasing sound speed with compression over a very limited range of compression ratios and only for unary phases, while the four binary phases all exhibi t increasing sound speed with compression down to V/V 0 = 0.78. Overall, the
89 Cheng potential exhibits better compression behavior than the Mendelev potential over the range 0.75 < V / V 0 < 1 Therefore, the Cheng potential is better suited for shock compressi on simulations and is used in the MD simulations presented in this work. 5.2.2 Amorphous Structure Creation T he compositions used to create the MG samples used in this work are Cu 30 Zr 70 Cu 50 Zr 50 and Cu 70 Zr 30 The initial fcc crystal structures consist of 4000 Cu atoms with 3D periodic boundary conditions. To achieve the desired composition, Cu atoms are randomly replaced with Zr atoms. Each system is given an initial temperature of 300 K and then heated at a constant rate of 10 K/ps to 2300 K. This temperature is chosen as it is far above the equilibrium melting temperatures of both Cu and Zr  After heating, the system is equilibrated for 200 ps at 2300 K to ensure uniform mixing. The system is then cooled to 300 K at a constant rate of 5 K/ps. The cooling rate is within an order of magnitude of those used in previous studies [91,93,211] The heating and quenching sim ulations are performed under the isobaric canonical ( NPT ) ensemble with a target pressure of zero. Mean square displacement ( MSD ) curves and radial distribution functions ( RDF s) are calculated for each composition after heating to 2300 K and after cooling to 300 K. The MSD curves of the system at 2300 K, shown in Fig. 5 2(a), exhibit linear behavior typical of a liquid. In contrast the MSD curves for the quenched samples, shown in Fig. 5 2(b), are four orders of magnitude smaller, which indicates the atoms are capable of only local rearrangements, as is true in solids. Upon cooling, the second peak of the RDF splits indicating the final solid structure is amorphous, see Fig. 5 2(c) and 5 2(d). The average atomic volume as a function of temperature V(T) during cooling is shown in Fig. 5 3 for all three compositions. Linear fits are applied to V(T) on either side of the
90 glass transiti on. The temperature at the intersection of the two fits gives an estimate of the glass transition temperature T g Figure 5 2. MSD s for (a) liquid and (b) metallic glass and Cu Cu RDF s for ( c) liquid and (d) metallic glass Table 5 1 gives the equilibrium density, elastic constants sound speeds and glass transition temperatures T g for each of the three MG samples. In the table, is the mass density, B is the bulk modulus, G is shear modulus, v is the Poisson ratio and c B and c L are the bulk an d longitudinal sound velocities, respectively. Systems with higher Cu content have higher values of T g which is consistent with experimental results  and previous MD simulations 
91 Figure 5 3. Volume tempe rature relationships during cooling for Cu 30 Zr 70 Cu 50 Zr 50 and Cu 70 Zr 30 systems Table 5 1 Calculated e lastic constants and glass transition temperature T g for Cu 30 Zr 70 Cu 50 Zr 50 and Cu 70 Zr 30 metallic glass samples. Composition (g/c c ) B (GPa) G (GPa) v c B (km/s) c L (km/s) T g (K) Cu 30 Zr 70 6.918 95.62 23.38 0.388 3.72 4.31 657 Cu 50 Zr 50 7.298 107.14 24.37 0.396 3.83 4.34 709 Cu 70 Zr 30 7.782 115.76 26.52 0.395 3.86 4.40 753 Note: densities, elastic constants and sound speeds are obtained for zero temperature. 5.2.3 Piston driven Shock Method To produce the large scale systems required for shock simulations, the small size MG samples are replicated along the x y and z directions. The + x axis is treated as the shock direction. The final cross sectional s iz es of the replicated systems are approximately 2020 nm 2 Cross sections of 4.24.2, 8.48.4 and 16.816.8 nm 2 were simulated in Ref.  and no size effects were found on the plastic deformation
92 behavior. The larger cros s section of our samples is mainly to provide better averages of thermodynamic properties extracted during shock compression simulations. The thickness of the sample along the x axis is 800 nm when U p when U p nm thick samples consist of ~20 million atoms, while the 400 nm thick samples consist of ~10 million atoms. Thicker samples are used for U p 1.5 km/s to allow for more time to track the evolution of the shock profile. At larger values of U p the shear stress goes to zero within a very s hort distance (< 50 nm) from the shock wave front, so thinner samples can be used. Large scale s amples are equilibrated within the NPT ensemble at 300 K and zero pressure for 50 ps to remove possible art i facts from the replication process. After equilibrat ion, the boundaries along the x axis are set to free that is there are free surfaces at the two ends; boundaries along the transverse dimensions remain periodic. Shock waves are generated in the equilibrated samples using the piston method  A thin slab (2 nm) located at the left edge of the system is treated as the piston. At the star t of the simulation, atoms inside the piston are assigned a uniform velocity along the + x axis and zero velocity along the transverse axes. Forces are zeroed for atoms within the piston at each time step. The piston velocity is ramped from 0 to U p within 1 0 ps, thereby generating a compression wave within the sample. The final piston velocity U p is varied from 0.125 to 2.5 km/s to obtain shocks of increasing strength. Shock simulations are carried out in the microcanonical ( NVE ) ensemble with a 1 fs time st ep. Material response during shock wave propagation is monitored using one dimensional profiles of thermodynamic properties. At the start of the simulation, the sample is divided along the x axis into small 1 nm thick slabs (~30,000 atoms). During
93 the simu lation, per atom thermodynamic quantities are summed within each slab and over multiple time steps. The sums are then appropriately normalized to get one dimensional profiles, which are outputted every 0.5 ps. Thermodynamic quantities used in the calculati on of the profiles include the pressure tensor P ij and temperature T where T is the kinetic temperature after subtracting out the center of mass velocity of the slab along the x axis. The shock pressure is P xx and shear stress is =0.5( P xx P yy )  where P yy P zz for the MG samples used in our simulations. 5.3 Sim ulation Results 5.3.1 Shock Wave Structure at Different Shock Intensities Profiles of shock pressure P xx and shear stress for Cu 50 Zr 50 MG at three different piston velocities U p are shown in Fig. 5 4. For U p = 0.125 km/s, there is no obvious drop in shear stre ss following initial shock compression and the shear stress remains almost constant behind the shock wave front. Because no appreciable plastic deformation occurs at this stress level, the shock wave can be considered elastic. For U p = 0.5 km/s, there is s ignificant relaxation of shear stress within the shock, which indicates the MG has undergone plastic deformation and the final pressure state is above the Hugoniot elastic limit ( HEL ) The shear stress reaches its peak value peak rapidly following the initial compression by the elastic precursor, then drops as plastic deformation starts to develop. This is in contrast to previous MD shock compression simulations of Cu 54 Zr 46 MG that used the Mendelev potential where the start of p lastic deformation did not necessarily coincide with the peak shear stress  Tracking the position of the peak shear stress, a demarcation between the elastic precursor and start of the plastic wave front in the shock pressure profile can be identified, as shown in Fig. 5 4. Behind the
94 plastic shock wave front, the shear stress continues to relax as plastic flow develops. The f low stress flow is defined as the average quasi steady shear stress behind the plastic wave front. Figure 5 4. Shock wave profiles of sho ck pressure and shear stress from MD simulations of Cu 50 Zr 50 metallic glass at U p = 0.125, 0.5 and 1.5 km/s. Previous MD simulations of shock compression of metallic single crystals near the HEL show clearly defined elastic and plastic wave structures with a clear separation between the elastic and plasti c shock wave fronts [83,85,203,204] However, for MGs it is difficult to distinguish the plastic wave front and elastic precursor clearly, as the plastic wave front is smooth and sluggishly rising with no clear separation from the
95 leading elastic precursor. Such behavior is similar to what has been observed in small grain size nanocrystalline metals where grain boundary motion dominates under shock compression [84,207,208,228] At a much higher piston velocity, U p = 1.5 km/s, the plastic shock wave has effectively overdriven the elastic precursor, which retains a small but finite width. Behind the shock wave front, the shear stress drops to nearly zero within a very short distance (~4 0 nm). 5.3.2 Hugoniot Elastic Limit At the HEL, the material transitions from a pure elastic state to an elastic plastic state. The shock pressure at the HEL is denoted by P HEL Equivalently, P HEL is the amplitude of the elastic precursor when both elastic and plastic shock waves are present. As was already stated, the peak shear stress coincides with the onset of plastic deformation in our simulations. Thus, the peak shear stress also coincides with the location of the amplitude of the elastic precursor. Using the x position of the peak shear stress, the corresponding amplitude of the elastic precursor can be obtained from the shock pressure profile as shown by the vertical dotted line in Fig. 5 4. This pressure P HEL is shown by the blue horizontal dotted line. As the shock wave propagates through the sample, there is a gradual decay in P HEL of ~510 3 GPa/nm. The observed decay rate is virtually independent of the MG composition. Figure 5 5 shows P HEL as a function of the final shock pressure for Cu 30 Zr 70 Cu 50 Z r 50 and Cu 70 Zr 30 MGs. In addition to the gradual decay in P HEL there are also small fluctuations of ~0.1 GPa. To account for these effects, the data plotted in Fig. 5 5 represent averages over several profiles for a simulation at a given final shock pressure (or U p ) taken within a short time window of several ps. All three compositions show a similar linear dependence of P HEL on the final shock pressure. In contrast, pr evious MD
96 shock compression simulations of Cu 46 Zr 54 MG that used the Mendelev potential showed a constant P HEL of ~7 GPa independent of the final shock pressure  The pressure independent HEL found for the Mendelev potential may be attributed to the unphysical behavior of the sound speed with compression shown in F ig. 5 1. Due to this behavior, it is possible the leading elastic wave does not steepen into a shock wave, which leads to the onset of plastic deformation at the same pressure level regardless of the final shock pressure. We define the minimum P HEL to be t he lowest shock pressure that results in plastic deformation behind the shock wave front subject to the length and time scale limitations of MD. Minimum P HEL values measured in this work are 6.4, 7.0 and 7.3 GPa for Cu 30 Zr 70 Cu 50 Zr 50 and Cu 70 Zr 30 MGs, re spectively. Figure 5 5. Shock pressu re at the HEL P HEL as a function of final shock pressure for Cu 30 Zr 70 Cu 50 Zr 50 and Cu 70 Zr 30 metallic glass 5.3.3 U p U s Shock Hugoniot The time evolution of the shock pressure profile in Cu 50 Zr 50 MG for two different piston velocities U p = 0.5 and 1.5 km/s is shown in Fig. 5 6 The initial compression
97 wave created in the MG sample steepens as it propagates, which can be seen clearly in the inset of Fig. 5 6 (a) as a converging of the x t contour lev els starting from t = 0. At ~25 ps, the contour levels have converged, which indicates that the compression wave has broken into a shock. Parts of the sample that have not yet been affected by the shock wave have a pressure near zero and are shown in blue in Fig. 5 6 while areas that have been compressed to the final shock pressure are shown in red. For the overdriven plastic shock wave shown in Fig. 5 6 (a), the boundary region between these two states is very narrow and can be approximated by a single lin e, the slope of which roughly corresponds to the shock velocity U s A similarly narrow boundary region is also observed for elastic shock waves with amplitudes < P HEL Figure 5 6 Shock pressure distributions at different times and positions for Cu 50 Zr 50 metallic glass at (a) U p = 1.5 km/s and (b) U p = 0.5 km/s. However, when there is an elastic plastic two wave structure, as is this case for U p = 0.5 km/s, there is ambiguity in the location the plastic shock wave front. This ambiguity can be attributed to the sluggishly rising nature of the plastic front, which is evident from the wide transition region between the unshocked and shocked regions of
98 the sample in Fig. 5 6 (b). By using an appropriate value of the shock pressure P xx the positions of the elastic and plastic fronts as a function of time are extracted from the simulation. The elastic and plastic shock wave velocities are then estimated from the slopes of the corresponding trajectories. To track the elastic precursor, a P xx level with a value < P HEL is used. For Cu 50 Zr 50 MG, P HEL is ~11 GPa when U p = 0.5 km/s. Taking P xx levels at 5.7 and 8.3 GPa gives elastic precursor velocities of 4.6 53 and 4.649 km/s, respectively, which is a difference of < 0.1%. Therefore, we conclude that the elastic precursor velocity is insensitive to the exact choice of the P xx tracking level provided it is < P HEL This insensitivity can be attributed to the sharp rise and narrow width (~5 1 0 nm) of the elastic precursor. In contrast, the plastic shock wave velocity is sensitive to the choice of the P xx level used to track its motion through the sample For example, taking P xx levels of 13.52 and 15.26 GPa for the plastic front, gives velocities of 4.54 and 4.32 km/s, respectively, which is a difference of 4.8%. To calculate the plastic shock wave velocity in a way that is consistent for different shock intensities and compositions, we set the P xx level to trac k the plastic front equal to the median pressure in the plastic shock ( P shock + P HEL )/2, where P shock is the final steady shock pressure and P HEL is the pressure at the peak shear stress Note P HEL depends on the final shock pressure and MG composition as is shown in Fig. 5 5. The U p U s relationships for shock compression of Cu 30 Zr 70 Cu 50 Zr 50 and Cu 70 Zr 30 MGs obtained from MD are plotted in Fig. 5 7 The U p U s plots shown in Fig. 5 7 can be broken down into three regions. When U p < 0.25 km/s (region I), shear stress relaxation behind the shock wave front is negligible and a single elastic shock wave develops with an amplitude < P HEL As U p is
99 increased beyond 0.25 km/s, plastic deformation occurs behind the elastic shock wave fron t and a plastic wave develops. Within the interval 0.25 < U p < 0.75 km/s (region II), the plastic wave velocity is less than that of the elastic precursor and the two waves separate during propagation. However, the rate of increase of the plastic wave velo city within this interval is greater than that of the elastic precursor. When U p > 0.75 km/s (region III), both the elastic precursor and plastic wave have the same velocity and propagate as a single elastic plastic shock wave. Figure 5 7 Shock velocit y U s as a function of piston velocity U p for Cu 30 Zr 70 Cu 50 Zr 50 and Cu 70 Zr 30 metallic glass The U p U s data for each MG composition are found to exhibit approximately piecewise linear behavior and are fit to relationships of the form U s = c 0 + sU p where c 0 and s are fitting parameters. Each elastic branch is best described by a single linear relationship, while each plastic branch is best described by two different linear
100 relationships: one for region II and another for region III. Table 5 2 gives the best fit parameters for both the elastic and plastic branches of the U p U s Hugoniot for Cu 30 Zr 70 Cu 50 Zr 50 and Cu 70 Zr 30 MGs. Both c 0 and s increase with increasing Cu content for all three regions. The increase in c 0 is expected as c 0 c B for the plastic bran ch and c 0 c L for the elastic branch and, as can be seen in Table 5 1 both c B and c L increase with increasing Cu content. Table 5 2 Calculated parameters f or both the elastic and plastic branches of the U p U s Hugoniot for Cu 30 Zr 70 Cu 50 Zr 50 and Cu 70 Zr 30 metallic glass Composition Elastic (regions I & II) Plastic (region II) Plastic (region III) c 0 (km/s) s c 0 (km/s) s c 0 (km/s) s Cu 30 Zr 70 4.26 0.42 3.60 1.4 2 3.9 7 0.85 Cu 50 Zr 50 4.3 5 0.62 3.8 1 1.4 3 4.1 3 0.90 Cu 70 Zr 30 4.4 3 0. 80 3.96 1.52 4.25 1.02 Note: t he standard errors in c 0 and s are < 1% and < 7%, respectively. 5.3.4 Composition Dependence of the Shock Response Shock wave profiles at U p = 1.5 km/s for Cu 30 Zr 70 Cu 50 Zr 50 and Cu 70 Zr 30 MG are shown in Fig. 5 8 Profiles are taken at 100 ps after the initial piston compression. As Cu content increases from 30 % to 70%, the shock pressure increases from 54 to 66 GPa. This is consistent with the increase in bulk modulus with increasing Cu content of the MG, as shown in Table 5 1 Figure 5 8 also shows tha t peak increases towards Cu rich MG compositions. The dependence of the stress state on Cu content is nonlinear. For example, at U p = 1.5 km/s, peak = 2.7 GPa for Cu 30 Zr 70 while peak = 3.2 GPa for
101 Cu 50 Zr 50 which is an increase of ~18%. However, peak = 4.7 GPa for Cu 70 Zr 30 which is an increase of 47% over that of Cu 50 Zr 50 Figure 5 8 Shock wave profiles at t = 100 ps for Cu 30 Zr 70 Cu 50 Zr 50 and Cu 70 Zr 30 metallic glass samples at U p = 1.5 km/s. The temperature of the material in the aftershock flow is plotted in Fig. 5 9 as a function of the final shock pressure. In the figure, the d ashed vertical line indicates departure from linear behavior at low shock pressures. For shock pressures < 25 GPa, the temperature increases almost linearly. This is true for all three MG compositions. The temperature increase from elastic deformation is independent of composition and is modest. At shock pressures > 25 GPa, the temperatur e starts to show significant deviations from linearity. As we know, the temperature reflects the intensity of the
102 movement of atoms and local structures (STZs for MGs). The higher the temperature, the more intense the activities of STZs for MGs. The shock induced temperature change decreases with increasing Cu content for the same shock pressure. When the shock pressure is around 60 GPa, the temperature of Cu 30 Zr 70 MG is ~130 K higher than Cu 70 Zr 30 MG. This indicates that Zr has a more important effe ct on t emperature rise than Cu. Figure 5 9 Shock induced temperatures a s a function of shock pressure. Figure 5 1 0 shows peak and flow as a function of the final shock pressure. The peak shear stress coincides with the start of plastic deformation in our s imulations provided the final shock pressure > P HEL Thus, the peak shear stress peak is a measure of the critical shear stress necessary to initiate plastic flow within the pristine glass under compress ion by the elastic precursor. The flow stress flow is the average of the residual shear stress far away from the plastic shock wave front a nd gives a measure of the strength of the material within the aftershock flow.
103 The peak shear stress increases with increasing shock pressure (states A C in Fig. 5 10 ( a)) until reaching a maximum value (state D) that is dependent upon the composition of the MG. Beyond this maximum, peak decreases indicating softening of the MG. As peak corresponds to the shear strength of the initial elastically compressed MG structur e without existing STZs to facilitate plastic flow, the softening must be thermal in nature. The higher the Cu content of the MG, the higher the maximum value of peak and the higher the shock pressure required before the MG starts to soften. For example, peak reaches a maximum value of ~2.8 GPa at a shock pressure of ~53 GPa for Cu 30 Zr 70 MG, whereas peak reaches a maximum value of ~4.8 GPa at a shock pressure of ~82 GPa for Cu 70 Zr 30 MG. Figure 5 10 (a) Peak shear stress and (b) flow stress as a funct ion of the final shock pressure The flow stress flow increases at low shock pressures (state A in Fig. 5 1 0 (b)), where the shock deformation is principally elastic, then reaches a maximum value at the onset of plastic deformation (state B). A t a given shock pressure flow increase s with
104 increasing Cu content, i.e., Cu 30 Zr 70 and Cu 50 Zr 50 MGs both yield more easily than Cu 70 Zr 30 As was discussed about Fig. 5 8 the dependence of the stress state on composition is nonlinear. For example, the maximum valu e of flow for Cu 30 Zr 70 and Cu 50 Zr 50 MGs are ~0.7 and ~0.8 GPa, respectively, while for Cu 70 Zr 30 flow has a maximum value of ~1.0 GPa. B eyond its maximum value, flow starts to decrease with increasing shock pressure (state C in Fig. 5 1 0 (b)). Such soften ing was also observed in previous MD simulations of shock compression of Cu 46 Zr 54  and Cu 46 Zr 46 Al 8  MGs, where it was attributed to a combination of thermal softening due to shock heating and strain softening due to formation of STZs and subseque nt strain localization. At high shock strengths, flow goes to zero (state D) indicating a complete loss of strength for the material within the aftershock flow. The point at which this transition occurs is sensitive to the MG composition. For Cu 30 Zr 70 Cu 50 Zr 50 and Cu 70 Zr 30 loss of flow strength occurs when the shock pressure (temperature) is ~67 GPa (1400 K), ~73 GPa (1500 K) and ~100 GPa (2200 K), respectively. Clearly, increased Cu content leads to an increased resistance to plastic flow that persists up to very high shock pressures. The Cu Cu RDFs for Cu 50 Zr 50 corresponding to states A D are shown in Fig. 5 11 The RDFs are calculated for a 50 nm slice of material along the x axis far behind the plastic shock wave front (not at the peak stress state). For states A C, the glassy structure is largely maintained after shock compression, while for state D the split in the second peak has been reduced significantly indicating a tra nsformation to a more liquid like structure. However, from this analysis alone, it is not clear whether the change in
105 the local structure can be attributed to melting. The issue of shock induced melting will be examined in detail in Chapter 6. Figure 5 1 1 Cu Cu RDFs in the aftershock flow of Cu 50 Zr 50 metallic glass for fo ur different shock intensities. 5.4 Discussion 5.4.1 Pressure dependent Yield Criterion Our MD simulations show, independent of Cu composition, that P HEL increases with increasing shock intensity see Fig. 5 5 The mean pressure P m = ( P xx +2 P yy )/3 at the HEL increases with increasing shock strength, therefore P HEL increases with the magnitude of the confining stress. Several impactor based shock experiments on multi component MGs have also observed a similar pressure dependent HEL [57,60,200] Figure 5 12 shows P HEL as a function of P m at the HEL for our MD simulations along with experimental data from Refs. [57,60,200] In experiment, transverse pressure info rmation is not available, so P m
106 error in P m of ~5%. Good agreement between our MD simulation data and shock compress ion experiments on a variety of different multi component MGs suggests that a pressure dependent HEL is a universal feature of MGs. Figure 5 1 2 Shock pressure at the HEL P HEL versus mean pressure P m at the HEL. Pressure dependent yield behavior has been observed in MGs under many different loading conditions [229 235] and several different variants of yield crite ria for MGs that account for the effect of non shear components of the pressure tensor have been proposed  Lu and Ravichandran found that the yield strength of Zr 41.2 Ti 13.8 Cu 12.5 Ni 10 Be 22.5 MG (Vitreloy 1  ) measured in quasistatic compression tests was well described by a pressure dependent Tresca criterion of the form Y = Y 0 + P m with = 0.34  In our MD simula tions, the yield strength Y = 2 peak where peak is the peak shear stress obtained from the shock profiles. Fitting the yield strength data as a function of P m at the yield point for our simulated MGs, we obtain = 0.26 0.3. The difference in between s imulation and experiment may be attributed to different
107 compositions of MGs as well as differences in sample preparation, specifically the fact that the cooling rate used to form the MGs in our MD simulations is many orders of magnitude higher than in expe riment. 5.4.2 Bilinear U p U s Hugoniot Our simulated U p U s Hugoniot data exhibit a plastic branch characterized by approximately bilinear behavior with a change in slope occurring at the point when the plastic wave overtakes the leading elastic precursor, i.e., at the transition from region II to region III in Fig. 5 8. The sh ock response of Zr 51 Ti 5 Ni 10 Cu 25 Al 9 MG was investigated in two separate experiments [65,200] Both experiments found an elastic plastic response with P HEL ranges that were consistent with one another. The plastic branches of the U p U s Hugoniot curves for both experiments were found to exhibit good linearity. However, when both sets of U p U s da ta are combined, a kink appears at a shock pressure of ~18 GPa ( U p ~ 0.55 km/s). In our MD simulations, the transition from region II to region III occurs at a shock pressure of ~ 25 GPa ( U p ~ 0.75 km/s). On the one hand, the fact that both experiments and simulation agree on the existence of a change in slope in the U p U s Hugoniot data suggests that our results may not necessarily be an art i fact of the limited length and time scales of MD simulation. On the other hand, previous MD simulation results for C u 46 Zr 54 MG showed no such change of slope with the plastic branch exhibiting good linearity throughout the entire range of shock pressures  Therefore, further work is needed to determine the origin of the change in slope observed in the plastic branch of our simulated U p U s Hugoniot data. 5.5 Summary The shock respons e of Cu x Zr 100 x ( x = 30, 50 and 70) MG has been investigated over a wide range of shock intensities using large scale MD simulations. A two wave
108 elastic plastic structure is found for all three compositions. At low shock pressures, the plastic wave propagates more slowly than the leading elastic shock wave and the two waves separate over time. As shock pressure is increased, the plastic wave eventu ally catches up to the elastic precursor and the two form a single shock wave structure. Within the split wave and overdriven regimes, the amplitude of the elastic precursor P HEL is found to increase with increasing confining stress, thereby indicating a p ressure dependent yield criterion for the simulated Cu x Zr 100 x MGs. The U p U s data for the plastic branch of the Hugoniot were found to exhibit approximately bilinear behavior with the slope gradually decreasing as the velocity of the plastic wave approach es that of the elastic precursor. The Cu content of the MG is found to have a significant effect on the shock induced elastic plastic transition with resistance to plastic deformation increasing towards Cu rich compositions. Specifically, MGs with higher Cu content can support larger shear stresses before yielding both within the pristine glass and within the aftershock flow. Previous MD simulations  have shown such composition dependence to be a result of changes to the local structure of the MG. In Chapter 6, the microscopic deformation behavior of Cu x Zr 100 x MGs under shock loading will be examined in detail. Emphasis will be placed on how such changes depend on the intensity of the shock and the Cu content of the MG.
109 CHAPTER 6 MICROSTRUCTU RAL RESPONSE OF METALLIC GLA SSES TO SHOCK LOADING 6.1 Background The Cu Zr system has served a prototypical role in both exp eriment an d MD simulations owing to its good glass forming ability over a wide composition range [216,218,237 239] and availability of interatomic potentials [217,240,241] Often new insights gained from studying this system are applicable to other metal met al binary MGs as well as Zr and Cu based multicomponent MGs. Both experiments and MD simulations have shown that the mechanical behavior of Cu Zr MGs is sensitive to composition. Quasistatic uniaxial compression experiments found increasing yield strength and decreasing plastic strain as Cu content increases from 46 % to 65%  Moreover, MD simulations of shearing, uniaxial and tensile loading of Cu Zr MGs at high strain rates found that the yield strength and the propensity for strain localization increase as Cu content increases from 30 % to 70%  The fundamental unit of plastic deformation in MGs is the shear transformation zone (STZ), wh ich is a small cluster of atoms that undergoes a shear event in response to an applied shear stress [2,47,242 245] Even in monolithic glasses, there are fluctuations and variations in the short range order (SRO) that lead to regions of l ow elastic stiffness or low stability [245 248] Such sites serve as prime candidates for STZ initiation. The structure of C u Zr MGs is best understood from the perspective of efficient packing of Cu centered clusters  MD simulations have shown that as Cu content increase s, the average coordination numbe r around Cu increases towards twelve and the fraction of Cu centered icosahedral clusters grows dramatically [217,219] Another MD study found that regions of high icosahedral order within
110 Cu 50 Zr 50 MG were more resistant to plastic flow than regions where icosahedral order was low  Shock waves have long been used as a tool to measure the dynamic strength of materials at high strain rates > 10 5 s 1 While several impactor based shock experiments have been performed on a variety of different MGs [57 64] measurements performed during shock propagation are limited to surface displacements or velocities, and microstructural data relate d to internal deformation of the sample can only be obtained after shock release. However, MD simulations have access to the complete dynamic response of the system at all stages of shock compression, while overlapping with the upper range of strain rates attainable in experiment. Therefore, they offer a unique opportunity to study the microscopic deformation behavior of MGs at high strain rates. A few MD simulations have already investigated the shock response of Cu Zr MGs [91,93,211] However, such studies focused on a single MG composition, so changes in the plastic deformation behavior due to changes in composition were not addressed. To this end we perform MD simulations to investigate the microstructural deformation behavior of Cu x Zr 100 x MGs under shock loading. This builds on the work done in Chapter 5, which focused mainly on the shock wave structure and Hugoniot states. We take x to be 30 % 50 % and 70% which corresponds to the same compositions used in Chapter 5 and also overlaps with the composition ranges used in the experiments  and MD simulations  mentioned above. We show that localized shear regions or STZs are the main carrier o f plastic deformation in the simulated MGs. Near the onset of elastic plastic behavior, STZ activity is inhomogeneous, but becomes more homogeneous as shock intensity increases. Shock
111 induced melting is shown to occur at the highest shock intensities. Plas tic strain within the aftershock flow is found to decrease with increasing Cu content of the MG. Changes in the SRO of the MGs are analyzed from the perspective of Cu centered clusters. Icosahedral clusters are shown to be the most stable during shock comp ression. 6.2 Simulation Methods MD simu lations are performed using LAMMPS  Interatomic interactions are described by an EAM potential for the Cu Zr( Al) system developed by Cheng et al.  In contrast, previous MD studies on shock compression of Cu Zr MGs [91,93] used the Finnis Sinclair EAM potential developed by Mendelev et al.  Following the work in Chapter 5, t he Cheng potential is chosen as it exhibits good behavior under compression. Initial structures consist of a 126.96.36.199 n m 3 cube of ~24,000 Cu and Zr atoms arranged in the B2 (CsCl type) structure with periodic boundary conditions along x y and z Depending on the desired composition, an appropriate fraction of Cu or Zr atoms is randomly replaced by atoms of the opposite type. Initial s tructures are heated at a constant rate of 10 K/ps t o a temperature of 2300 K. As was shown in Chapter 5, this temperature is high enough to ensure that the initial structures have completely melted After heating, the melts are equilibrated at 2300 K for 200 ps. Samples are then cooled at a rate of 10 K/ps to a final temperature of 300 K. The high cooling rate inhibits crystallization resulting in glass formation. The amorphous structure of the samples was con firmed using similar analysis to what was performed in Chapter 5. T he cooling rate of 10 K/ps is comparable to previous MD simulations [91,93,209,211] Glass formation simulatio ns are carried out in the isobaric isothermal ( NPT ) ensemble with a target pressure of zero and a time step of 2 fs.
112 To generate the large scale samples required for shock compression simulations, the amorphous structures obtained following cooling are rep licated along the x y and z directions. The shock direction is chosen as the + x axis. Sample dimensions after replication are 797.219.518.2, 798.222.821.9 and 799.122.621.4 nm 3 for Cu 30 Zr 70 Cu 50 Zr 50 and Cu 70 Zr 30 MGs, respectively. The large size al ong the x axis is used to allow more time for relaxation processes occurring behind the shock wave front to develop As was stated in Chapter 5, the cross sections of our samples are larger than those used in previous MD shock simulations [90,93] The as replicated samples are heated at 10 K/ps to 700 K and then annealed for 200 ps in the NPT ensemble at a target pressure of zero. The annealing step is performed to reduce the long range order in the sample produced by the replication process. After annealing, the samples are cooled at a rate of 10 K/ps to 300 K. Shock compression simulations are performed using the piston method described in Chapter 5. Material response during shock wave propagation is monitored using the 1D binning analysis described in Chapter 5 Profiles of the pressure tensor P ij and temperature T are calculated, where T is the kinetic temperature subtracting out the center of mass velocity of each slab along the x axis Shear stress = 0.5( P xx P yy ) profiles are also calculated, where P yy P zz for the MG samples used in our simulations. The mean square displacement (MSD) of atoms in the aftershock flow is calculated to test for shock induced melting. Because of the nonzero flow velocity of atoms along the x axis, only displacements along the y and z dimensions are used in the calculation of the MSD.
113 Plastic d eformation of the MG during shock compression is monitored using the von Mises shear strain i vM  Calculation of i vM requires two atomic configurations, one current, and one reference. The sample after annealing is used as the reference atomic configuration. Variations of < 1 GPa and < 3% are found in the local pressure tensor and composition for the reference configurations for each of the three MG compositions. To calculate i vM f irs t a local transformation matrix J i is determined that best maps where is the separation (row) vector between atoms i and j ( 0 denotes the reference configuration) In Eq. 6 1, atom j is one of atom i nearest neighbors and is the total number of nearest neighbors of atom i at the reference configuration within some cutoff distance In all our calculations, the cutoff distance is chosen as 4.5 . The local transformation matrix J i for each atom is determined by minimizing The strain matrix i for each atom is computed as where T represents the transpose and I is the identity matrix. The local von Mises shear strain i vM for each atom is then given by ( 6 1 ) ( 6 2 ) ( 6 3 ) ( 6 4 )
114 where the index i in the components of the strain matrix i is omitted for clarity. The SRO of the MG is characterized using a polydisperse Voronoi tesse llation analysis  wh ere the atomic radii for Cu and Zr are 128 and 160 pm, respectively. The Voronoi polyhedron associated with each atom is labeled by a set of indices < n 3 n 4 n 5 n 6 >, where n i is the number of faces with i vertices. For example, the Voronoi polyhedron with indices <0,0,12,0> is composed of 12 faces each face having 5 vertices/edges. In a perfect fcc or bcc crystal, there is only a single unique Voronoi polyhedron with indices of <0,12,0,0> or <0,6,0,8>, respectively. When performing the Voronoi tessellation, any f aces having areas less than 1% of the t otal polyhedron surface area are deleted to reduce effects of degeneracy and thermal vibration  The coordination number CN of each atom is equal to the sum of the faces of its corresponding Voronoi polyhedron, i.e., V isualization of th e simulation snapshots, Voronoi tessellation analysis and calculation of i vM are performed using the OVITO software package  6.3 Simulation Results 6.3.1 Shock induced Deformation Shock compression simulations are performed for Cu 50 Zr 50 MG varying U p from 0.125 to 2.0 km/ s Figure 6 1 shows shock pressure P xx and shear stress profiles at t = 125 ps for U p = 0.125, 0.375, 0.875 and 1.75 km/s. At the lowest piston velocity, U p = 0.125 km/s, the shear stress level behind the shock wave front is nearly constant. For example, th e average shear stress in a 100 nm thick slab immediately behind the shock wave front is 0.50 GPa, while 550 nm behind the shock wave front the average shear stress is 0.47 GPa. Therefore, shear stress relaxation is minimal, and the shock can be considered elastic at this low piston velocity
115 Figure 6 1 Pressure P xx and shear stress profiles at t = 125 ps for Cu 50 Zr 50 metallic glass at different piston velocities. When U p is increased to 0.375 km/s, shear stress relaxation behind the shock wave front is apparent, which indicat es the material has undergone some plastic deformation. This is consistent with the results presented in Chapter 5 that showed the onset of plasticit y occurs near U p = 0.25 km/s. However, the relaxation process is slow and a significant residual shear stress, or flow stress flow remains far behind the shock front. As U p is increased further to 0.875 km/s, the shear stress relaxation becomes more rapi d and more complete, as indicated by the sharp drop in shear stress behind
116 the shock front and a smaller flow level. This softening trend continues with increasing shock strength. At U p = 1.75 km/s, the shear stress goes to zero within ~50 nm behind the s hock front. Figure 6 2 Snapshots of von Mises shear strain for Cu 50 Zr 50 metallic glass at different piston velocities. Position of the shock front in each snapshot is denoted by a vertical line. Snapshots of the atomic shear strain i vM for Cu 50 Zr 50 are shown in Fig. 6 2 for the same values of U p plotted in Fig. 6 1. In each snapshot, the location of the shock front is denoted by a vertical line. For the shock wave generated at U p = 0.125 km/s,
117 ma terial behind the shock front has no perceptible change in i vM relative to the ambient state ahead of the shock front. This confirms that the deformation is almost completely elastic. The average i vM in the ambient unshocked material is ~0.04, which is i n agreement with previous MD shock simulations in Cu 46 Zr 54 MG performed using the Mendelev potential  At U p = 0.375 km/s, a small number of STZs have nucleated behind the shock front. These STZs are isolated from one another and are embedded in a matrix of material with low i vM (near the ambient value). Further behind the shock front, there is some growth and interconnection of STZs, but, for the most part, they remain spatially separated. The average temperature s of the STZs are similar to the temperature of the surrounding mat rix, which indicates STZ nucleation is not due to thermal fluctuations. Instead, variations in the SRO of the MG dictate the observed STZ nucleation sites. This conclusion is in agreement with MD simulations of shock [91,93] and non shock [219,249,254 256] deformation of Cu Zr MGs. Figure 6 3 shows a snapshot of i vM for a 30 nm thick chunk of the sample containing the shock wave front along with profiles of vM ( x ) and ( x ). As can be seen clearly in the figure, the peak shear stress peak coincides closely with the onset of STZ nucleation. The part of the shock front from the beginning of the shock wave rise to the peak shear stress constitutes the elastic precursor. The start of the elastic precursor is visible as a small dip in vM ( x ), which is indicated by an arrow in the figure. At U p = 0.875 km/s, peak is nearly twice that at U p = 0.375 km/s. The high shear str ess level results in an avalanche of STZ nucleation events at the plastic front, which is consistent with the sharp relaxation observed in the shear stress profile shown in Fig.
118 6 1. Approximately 15 nm from the start of the plastic front the shear st ress has been reduced to half of its peak value. With increasing distance behind the plastic front, and thus time for evolution, the STZs grow larger and become more interconnected forming a dense network of plastic deformation. However, the deformation is not completely homogeneous and p ockets of low i vM (near the ambient value) are visible, which have larg ely resisted shear deformation. Figure 6 3 Snapshot at t = 40 ps of von Mises shear strain i vM and profiles of vM ( x ) and shear stress for Cu 50 Zr 50 and U p = 0.375 km/s. The i vM snapshot at U p = 1.75 km/s shows th at the shear strain distribution is now almost completely homogeneous behind the shock front. As was shown in Fig. 6 1, the shear stress drops rapidly from its peak value to zero within ~50 nm. The temperature immediately behind the shock front is ~1350 K and increases steadily to a
119 maximum of ~1550 K near the piston. As will be shown later, the MG has melted under the a ction of such a strong shock wave 6.3.2 Composition Dependence and Melting We now turn our focus to the composition dependence of the shock response. To aid in the analysis presented in the remainder of this chapter we find it useful to differentiate between t wo states within the shock: state A at the position of the peak shear stress and state B in the quasi steady flow far behind the start of the shock wave rise The she ar stress at states A and B are peak and flow respectively. Figure 6 4(a) shows peak a s a function of the mean pressure P m = ( P xx + 2 P yy )/3 at the peak shear stress, while Fig. 6 4(b) shows flow as a function of the final shock pressure. Figure 6 4. (a) Peak shear stress peak as a function of the mean pressure P m and (b) flow stress flow as a function of the final shock pressure P xx As was mentioned in Chapter 5, peak corresponds to the critical shear stress required to plastically deform the pristine MG an d is related to the yield strength Y = 2 peak Our simulations show that i nd ependent of composition peak increases with
120 increasing P m Such behavior is consistent with the results presented in Chapter 5. Interestingly, when P m peak depends only on P m and not on the composition of the MG. However, when P m > 10 GPa, the curves start to diverge due to more pronounced thermal softening with decreasing Cu content. After the initial STZ nucleation events at state A, the shear stress continues to relax as plasticity develops behind the plastic front. At state B, the shear stress has reached the quasi steady value flow which corresponds to the shear strength of the r ejuvenated or plastically deformed, MG. In our simulations, flow is calculated by averaging the shear stress in a 50 nm thick slab near the piston face over several ps when the shock is close to the opposite edge of the sample. At a given shock pressure, flow is greatest for Cu 70 Zr 30 followed by Cu 50 Zr 50 then Cu 30 Zr 70 Therefore, MGs with higher Cu content exhibit greater resistance to plastic flow in the rejuvenated state. This is consistent with the results of Chapter 5 and the behavior of Cu x Zr 100 x M Gs under pure shear loading found in previous MD simulations  For all co mpositions, flow decreases with increasing shock pressure indicating shear softening and thermal softening due to heating by the shock wave At high enough shock pressures (50 80 GPa), the MG completely loses its ability to resist shear, i.e., flow This transition occurs at ~53 GPa, ~65 GPa and ~80 GPa for Cu 30 Zr 70 Cu 50 Zr 50 and Cu 70 Zr 30 MGs, respectively. The differences in these values compared to those given in Chapter 5 may be due to differences in preparation of the MG samples, including quench rate, inclusion of the annealing step and size of the initial amorphous sample prior to replication. However, the ordering observed among the different compositions remains the same.
121 Figure 6 5. Snapshots of von Mises shear strain for Cu 30 Zr 70 Cu 50 Zr 50 and Cu 70 Zr 30 metallic glasses at P xx = 20 GPa. Snapshots of i vM at state B for P xx = 20 GPa are shown in Fig. 6 5 for all three MG compositions. Visual inspection of the snapshots shows that the amount of shear deformation is highest for Cu 30 Zr 70 follow ed by Cu 50 Zr 50 then Cu 70 Zr 30 Defining atoms having i vM > 0. 2 as participating in a shear transformation event  the fraction of atoms involved in STZs in Cu 30 Zr 70 Cu 50 Zr 50 and Cu 70 Zr 30 are 29.3%, 17.5% and 8.6%, respectively. The average i vM at state B is plotted in Fig. 6 6 as a function of shock pressure. Averages are taken over the same 50 nm thick slab used in the calculation of
122 Fig. 6 4(b). At all shock pressures, the average i vM is highest for Cu 30 Zr 70 and decreases with increasing Cu content. This is consistent with the composition dependent behavior of flow observed in Fig. 6 4(b) and confirms that composition plays a large role in the yield behavior of the MG Specifically, Cu x Zr 100 x MGs with higher Cu compositions (over the range of 30% t o 70%) exhibit greater resistance to shear deformation and less plastic strain in the aftershock flow. Figure 6 6. Average von Mises shear strain in the quasi steady flow state B behind the shock front as a function of the final shock pressure P xx The diffusivity D is calculated at state B to determine the melting transition under shock loading for each MG sample. For a given simulation, the 2D MSD for all atoms within a 50 nm thick slab near the piston is calculated over an interval of 80 ps when t he shock wave front is near the far side of the sample A linear fit is applied to the MSD( t ) data and D is proportional to the slope of the fit Figure 6 7 shows D as a function of the shock pressure P xx for our simulated MGs together with results from pr evious MD Hugoniostat simulations for Cu 46 Zr 54 MG obtained using the Mendelev potential  In
123 that work, the transition to shock melting occurs when P xx 60 GPa, which is evident by a rapid increase in D Our results for Cu 50 Zr 50 MG show a similar increase in D near ~ 60 GPa, thereby indicating a similar melting transition unde r shock loading for the Cheng and Mendelev potentials Figure 6 7. Diffusivity as a function of shock pressure P xx for Cu Zr metallic glasse s from this work and Cu 46 Zr 54 metallic glass from Hugoniostat simulations of Ref.  Similar increases in the diffusivity data are observed for Cu 30 Zr 70 and Cu 70 Zr 30 but at different sh ock pressures. A sharp increase in D occurs near ~50 GPa for Cu 30 Zr 70 whereas for Cu 70 Zr 30 the increase in D is more gradual and starts near ~70 GPa Clearly, the shock pressure required to induce melting increases as the Cu content of the MG increases over the range of 30% to 70%. The diffusivity analysis is consistent with the flow stress data plotted in Fig. 6 4(b).
124 6.3.3 Local Structural Changes During Shock Compression Previous experiments  and MD simulations [217,219] have found that the local structures of Cu Zr MGs are best understood from the perspective of Cu centered clusters. Each cluster can be represented by a Voronoi polyhedron with indices < n 3 n 4 n 5 n 6 >. Although faces with more than six vertices are observed in our samples, the fraction of polyhedra with more than six sided faces is negligible and all observed polyhedra can be distinguished using indices up to n 6 There are several hundred unique Cu centered polyhedra types that make up the initial MG samples. However, the majority belong to a small number of polyhedra types. Specifica lly, we identified 24 types that represent ~70% of the total po pulation for all t hree compositions. Out of those 24 polyhedra types, nine were selected based on the criterion of greatest population change before and after shock loading. The polyhedra types (from 1 to 9), listed in order of increasing CN, are <0,2,8,0>, <0,1,8,1>, <0,3, 6,1>, <0,2,8,1>, <0,1,8,2>, <0,0,12,0>, <0,2,8,2>, <0,1,10,2> and <0,3,6,4>. Together they make up ~40% of the overall polyhedra population of the initial MGs. Among these polyhedra are the so called Z clusters, which have the minimum intrinsic disclinatio n (rotational defect) density for a given CN and are presumably the most stable  These are <0,2,8,0>, <0,2,8,1>, <0,0,12,0> and <0,1,10,2> for CN = 10, 11, 12 and 13, respectively. Clusters with Voronoi indices <0,0,12,0> have no intrinsic disclinations and have the maximum possible five fold bonds for CN = 12. They are known as full icosahedra [259 261] Representative Cu centered clusters for each of the Z clust ers are shown in Fig. 6 8 in order of increasing CN from left to right. Cu atoms are colored brown, while Zr atoms are colored green.
125 Figure 6 8. Representative Cu centered clusters corresponding to the given set of Voronoi indices. Cu atoms are colored brown, while Zr atoms are colored green. We tracked the populations of the four Z clusters in a 50 nm thick slab near the piston as the shoc k wave passed through the sample Every 10 ps the populations we re calculated and compared to the values at the pre vious iteration When t populations reach ed nearly steady values. Figure 6 9 shows the fraction s of the different Cu centered Z clusters at their quasi steady values as a function of the p iston velocity U p for all three MGs. Independent of composition, the fraction of <0,2,8,0> polyhedra decreases with increasing U p while the fraction of <0,1,10,2> polyhedra increases up until U p This is consistent with the higher compression ratios at higher shock intensities, w hich leads to an increase in the average CN of the polyhedra. The behavior is somewhat more complicat ed for <0,2,8,1> and <0,0,12,0> polyhedra. Both Cu 30 Zr 70 and Cu 50 Zr 50 exhibit large increases in <0,0,12,0> polyhedra with increasing U p while in Cu 70 Zr 30 the fraction remains nearly steady. When U p the fraction of <0,0,12,0> polyhedra decreases due to shock induced melting.
126 Figure 6 9 Fraction of Cu centered Voronoi polyhedra at the quasi steady flow state behind the shock front as a functio n of the piston velocity U p The polyhedron for each coordination number (CN) corresponds to the Z cluster, which has the minimum disclinations of any polyhedra type for that CN. The data in Fig. 6 9 reflects total changes in polyhedra populations. Howeve r, it does not convey information related to transformations of polyhedra before and after shock loading. Such information would help assess the relative stability of different polyhedra types. To this end, we determined the distribution of polyhedra for a fixed group of atoms at two configurations: an initial configuration corresponding to a 50 nm thick slab in the unshocked sample near the piston and a final configuration corresponding to the same slab at the quasi steady flow state B behind the shock fro nt.
127 Voronoi tessellation was performed on both the initial and final atomic configurations. The fraction of polyhedra transformed into different types was then calculated by comparing the initial and final Voronoi indices for each Cu atom. For example, the fraction of <0,0,12,0> polyhedra in the initial state for Cu 50 Zr 50 MG is 5.7%. At U p = 0.75 km/s, 1.4% are still the same type at the final state, state B. Therefore, 0.014/0.057 = 24% of the original <0,0,12,0> polyhedra remain the same type, while 1 0.24 = 76% have transformed into other types. The fractions of transformed polyhedra a s a function of the piston velocity U p are plotted in Fig. 6 10 for the following nine polyhedra types: <0,2,8,0>, <0,1,8,1>, <0,3,6,1>, <0,2,8,1>, <0,1,8,2>, <0,0,12,0>, <0,2 ,8,2>, <0,1,10,2> and <0,3,6,4>. These polyhedra types were selected based on the criteria of highest overall population fractions and greatest total population change during shock compression. Data for each MG composition is shown. The polyhedra types in the legend are listed in order of increasing CN from top to bottom. Th e transform ed fraction can range from 100%, in which all the original polyhedra of that type have transformed into other types, to 0%, in which none have changed types. Note that the data in Fig. 6 10 only considers the initial and final atomic configurations, so mul tiple transformations may be occurring between these two end states. For all three compositions, the <0,0,12,0> polyhedra, or full icosahedra, exhibit the lowest transformed fraction over almost the entire range of U p This indicates that the full icosahed ra are the most resistant to change under shock loading. Such behavior is consistent with previous MD simulations of shock  and non shock [ 219] loading of Cu Zr MGs. The remaining eight polyhedra types, including the other three Z clusters,
128 exhibit less stability than the full icosahedra. Interestingly, when comparing the stability of the Z clusters: <0,2,8,0>, <0,2,8,1>, <0,0,12,0> and <0 ,1,10,2> with other clusters with the same CN, we find that the Z clusters are generally more stable. Figure 6 10 Fraction of Cu centered Voronoi polyhedra that have been transformed into other types after transitioning from the unshocked state to the q uasi steady flow state B behind the shock front. A value of 10 0% indicates that all initial polyhedra of that type have transformed into other types at state B, while a value of 0% indicates that none have changed types.
129 With increasing U p the fractions of transformed polyhedra increase then reach nearly stable values when U p increased more and more Cu centered clusters become involved in the plastic deformation processes behind the shock front. Then, when U p of transformed polyhedra saturate, corre sponding to an exhaustion of clusters participating in the plastic response Most polyhedra types shown in Fig. 6 10 exhibit such behavior apart from <0,0,12,0> polyhedra for Cu 70 Zr 30 MG, which continues to increase above U p we observe a gradual change in slope in the plastic branches of the U p U s Hugoniot curves between regions II and III that were shown in Chapter 5 6.4 Summary In this chapter we studied the microstructural changes that occurred during shock compression of Cu x Zr 100 x MGs using MD simulations. The main plastic deformation mechanism found in our simulations was n ucleation and growth of localized regions of higher than average shear strain, or simply STZs At low shock intensities, the number of STZ s nucleat ed behind the shock front is low and the plastic deformation is inhomogeneous. However, as the shock intensit y is increased, more STZs are nucleate d resulting in a more homogeneous plastic response MGs with higher Cu content exhibit higher shear strengths and less shear strain in the quasi steady flow behind the shock front, which is consistent with the increase in yield strength and decrease in plastic strain found in non shock MD simulations  and quasistatic loading experiments  The shock induced melting transition is also found to occur at higher shock pressures as the Cu content of the MG is increased. Full icosahedra are
130 found to be the most stable under shock compression. The exhaustion of the number of Cu centered clusters that are transformed under shock loading may contribute to the gradual change in slope observed in the plastic branch of the U p U s relationship shown in Chapter 5.
131 CHAPTER 7 FERROELECTRICTY IN LANTHANUM BOROGERMANATE 3 7.1 Background LaBGeO 5 is a prototypical and well s tudied material from the borosilicate stillwellite family of compounds, RB(Si,Ge)O 5 technological interest in this material stems from its ferroelectric properties [98 102] and non linear optical properties [108,262,263] Moreover, LaBGeO 5 has good glass forming ability and it is possible to nucleate nano crystallites within LaBGeO 5 glass  Such a ferroelectric nanocomposite material is a transparent and optically active medium and has interesting potential applications. The measured values for the spontaneous polarization are some what varied depending on the measurement method: pyroelectric measurements give a room temperature polarization of 4.2 C/cm 2  whereas calorimetric measurements give a saturated polarization of 2.7 C/cm 2  and second harmoni c generation intensity measurements give a value of 12 C/cm 2  At about 530C, LaBGeO 5 undergoes a second order transition to a high symmetry paraelectr ic phase [98 101] Valence Force Field examination of the phase transition  identified a possible phonon mode driving the ferroelectric phase transition. However, from an ab initio perspective, only examination of the low temperature (LT) paraelectric phase has been performed  ; neither the high temperature (HT) phase nor the dynamics of the phase transition have been examined. Consideration of both phases is necessary to determine the ferroelectric properties of LaBGeO 5 specifically the spontan eous polarization. In this chapter we pres ent DFT 3 The work described in this chapter has been published in B. J. Demaske, A. Chernatynskiy, and S. R. Phillpot, J. Phys: Condens. Matt. 28, 165901 (2016 ).
132 calculations on both the LT and HT phases of LaBGeO 5 as well as on the dynamics of the phase transition and the spontaneous polarization. We compare our results to available experimental data and previous theoretical models. 7.2 Simulation Methods The energetics and ferroelectric properties LaBGeO 5 a re determined from DFT calculations [126,265] Specifically, we use the pl ane wave DFT code VASP [131,159] which employs the PAW method to simplify treatment of core electrons [160,161] In our calculations, 11 electrons for La (5 s 2 5 p 6 5 d 1 6 s 2 ), 3 for B (2 s 2 2 p 1 ), 14 for Ge (3 d 10 4 s 2 4 p 2 ) and 6 for O (2 s 2 2 p 4 ) a r e explicitly treated as valence, while the remaining core electrons a re treated using the PAW pseudo potentials supplied by VASP For comparison purposes, the Perdew Zunger parameterization of the LDA  the PBE form of the GGA [12 9] and a revi sed PBE functional for solids PBEsol  a re used to approximate the exchange correlation functional. A plane wave energy cutoff of 550 eV centered Monkhorst Pack  333 k point mesh a re found to give energies converged to within 5 meV per atom. Us e of the LDA and GGA functionals results in an unphysical lowering of energy of the f orbitals of La such that they lie near the bottom of the conduction band. To correct for this behavior a Hubbard U correction i s applied to these orbitals using the Dudar ev formulation  present in VASP As there are no published U parameters avai lable for La in LaBGeO 5 we ch o ose instead to use published values of U = 11 eV and J = 0.68 eV for the f orbitals of La in LaTiO 3  This choice is motivated by the similar La bonding environments in La TiO 3 and LaBGeO 5 Structures a re relaxed until the maximum forces on each atom a re less than 0.01 eV/ and total stress on the unit cell i volume.
133 Second order force constants a re determined using the finite displa cement method implemented in the software package Phonopy  with VASP as the force calcula tor. Phonon dispersions ar e obtained by diagonalizing the dynamical matrix constructed from the second order force cons tants. The MEP between the two symmetry equivalent ferroelectric structur es i s calculated using the CI NEB method [168,169,192] with a total of seven intermediate images. The macroscopic polarization of LaBGeO 5 i s calculated using the modern theory of polarizati on based on the Berry phase formalism [119 123] as implemented in VASP. Figure 7 1 Structure of the low temperature phase of LaBGeO 5 from the side (left panel) and from the top (right panel). 7.3 Simulation Results 7.3.1 Low temperature Phase of LaB GeO 5 The LT phase of LaBGeO 5 belongs to the polar space group P 3 1 and contains three formula units per cell (24 atoms). A schematic of the crystal structure of the LT
134 phase showing the locations of the five unique O atoms is shown in Fig. 7 1. Each of the B and Ge ions are surrounded by oxygen forming BO 4 (blue) and GeO 4 (green) tetrahedra. La atoms are shown in gray and O atoms in red with indices labeled. The BO 4 tetrahedra are linked together through corners to form chains along the c axis, with each of the GeO 4 tetrahedra sharing corners with two BO 4 tetrahedra along a chain. The La ions are surrounded by oxygen forming large LaO 9 polyhedra that serve to interconnect the BO 4 /GeO 4 chains. These LaO 9 polyhedra are composed of 14 triangular faces with the O atoms arranged at the nine vertices as shown in Fig. 7 2. Table 7 1 Comparison of equilibrium lattice constants for the low temperature phase of LaBGeO 5 Experiment (20C) VASP (LDA+ U ) VASP (PBEsol+ U ) ABINIT (norm conserving) a () 7.0018 6.9986 ( 0.05%) 7.0563 (+0.78%) 6.957 ( 0.64%) c () 6.8606 6.8057 ( 0.80%) 6.8578 ( 0.04%) 6.742 ( 1.73%) c/a 0.9798 0.9724 0.9719 0.9691 Note: values in parentheses are percent age differences from experiment. Theoretical and experimental lattice constants for the low temperature phase of LaBGeO 5 are given in Table 7 1 Experimental values are taken from neutron powder data  Theoretical results from another DFT study using norm conserving pseudopotentials  are also given in addition to our DFT results P ercent age errors from experiment are given in parentheses. DFT results from our VASP calculations are given for the LDA and PBEsol function als with the Hubbard U correction. Results from the PBE functional are omitted as they overestimate experiment al lattice parameters; addition of the Hubbard U correction leads to an increase in the equilibrium volume,
135 worsening the agreement. The magnitudes of the total errors of the LDA+ U and PBEsol+ U lattice constants are almost equal. However, when considering the lattice constants from the high temperature phase as well as atomic positions and interatomic distances, PBEsol+ U outperforms LDA+ U in terms of overall accuracy. Therefore, the results presented in the remainder of this work a re calculated using PBEsol+ U functional unless stated otherwise. The DFT relaxed fractional coordinates of the atoms in the LT phase are given in Table 7 2 along with data from neutron powder experiments at 20C  and another DFT work  where norm conserving pseudopotentials were used in contrast to the PAW approach used here. Accounting for the differences in relaxed volumes, our atom positions differ from the neutron diffraction data by an average of < 0.1 with t he largest deviation of 0.08 for the O5 atom. Furthermore, our atom positions and those obtained using norm conserving pseudo potentials differ by on average ~0.1 with the largest difference of 0.37 for the B atom. Interatomic distances for the BO 4 an d GeO 4 tetrahedra and LaO 9 polyhedra are given in Table 7 3 F igure 7 2 Structure of the LaO 9 polyhedra (a) with attached GeO 4 (green) and BO 4 (blue) tetrahedra and (b) without.
136 T able 7 2 Comparison of atom positions for the low temperature phase of LaBGeO 5 Experiment (20C) This work (PBEsol+ U ) DFT (norm conserving) Atom Position and multiplicity x/a y/b z/c La 3a 0.4108 0.9990 1/3 0.0004 0.0010 0.0052 0.0019 0.0001 0.0013 B 3a 0.8847 0.9878 0.3086 0.0008 0.0035 0.0015 0.0460 0.0021 0.0017 Ge 3a 0.4181 0.0004 0.8433 0.0020 0.0068 0.0041 0.0010 0.0063 0.0071 O1 3a 0.1528 0.3428 0.0179 0.0024 0.0056 0.0027 0.0017 0.0046 0.0000 O2 3a 0.3356 0.1481 0.9989 0.0110 0.0069 0.0016 0.0081 0.0067 0.0065 O3 3a 0.1342 0.6089 0.3376 0.0031 0.0008 0.0085 0.0005 0.0007 0.0009 O4 3a 0.6081 0.1620 0.6615 0.0003 0.0029 0.0012 0.0018 0.0098 0.0051 O5 3a 0.0121 0.0543 0.7861 0.0070 0.0073 0.0089 0.0032 0.0036 0.0054 Note: theoretical values given as differences from experiment.
137 T able 7 3 Comparison of B O, Ge O and La O interatomic distances (in ) for the low temperature phase of LaBGeO 5 Experiment (20C) This work ( PBEsol + U ) DFT (norm conserving) La O () O1: 2.651 2.772 O1: 0.006 0.049 O1: 0.068 0.166 O2: 2.684 2.678 O2: 0.033 0.071 O2: 0.008 0.021 O3: 2.433 2.566 O3: 0.006 0.017 O3: 0.018 0.036 O4: 2.589 2.412 O4: 0.037 0.032 O4: 0.175 0.118 O5: 2.672 O5: 0.028 O5: 0.078 B O () O1: 1.480 O5: 1.414 O1: 0.027 O5: 0.042 O1: 0.021 O5: 0.034 O2: 1.535 1.481 O2: 0.025 0.021 O2: 0.029 0.035 Ge O () O1: 1.782 O3: 1.625 O1: 0.008 O3: 0.099 O1: 0.052 O3: 0.057 O2: 1.772 O4: 1.761 O2: 0.019 O4: 0.032 O2: 0.032 O4: 0.075 Note: theoretical values given as differences from experiment.
138 Figure 7 3. Electronic band structure and DOS for the low temperature phase of LaBGeO 5 calculated using DFT. The electronic band s tructure and DOS for the LT phase of LaBGeO 5 are shown in Fig. 7 3 In the figure, the z ero energy is shifted to the top of the valenc e band. Our DFT calculations predict LaBGeO 5 to have an indirect band gap of 4.06 eV between the without the Hubbard U correction, the calculated band gap is 4.12 eV. This value is in good agreement with the value of 4.54 eV determined from a previous DFT study  To the best of our knowledge, there are no e xperimental results for the band gap in LaBGeO 5 As was found in Ref.  M K concentrated group of bands near the upper edge of the conduction band belongs to the f orbitals of La. Without the Hubbard U correction, these bands are incorrectly predicted to lie near the bottom edge of the conduction band.
139 T he phonon dispersion and DOS for the LT phase of LaBGeO 5 is shown in Fig. 7 4. Negative numbers correspond to imaginary frequencies. The absence of any imaginary phonon modes indicates that the structure is mechanically stable. This agrees with experiment as there is no evidence of any further phase transitions in LaBGeO 5 for temperatures below 20C. Figure 7 4. Phonon dispersion and DOS for the low temperature phase of LaBGeO 5 calculated u sing DFT. 7.3.2 High temperature Phase of LaBGeO 5 At 532C, LaBGeO 5 undergoes a second order phase transition from the polar space group P 3 1 to the non polar space group P 3 1 21 [116,118] According to the model in Ref. [ 116] both the La and Ge atoms undergo small displacements to bring them to lie on the twofold axis, while the four unique oxygen (O1, O2, O3 and O4) in the LT phase merge into only two unique oxygen (O1 and O3) in the HT phase. The B atom
140 and the remai ning O5 atom lie near, but not on, the twofold axis. To maintain the stoichiometry of the cell, the generated sites on the B and O5 sublattices must be half occupied. It is not clear whether the disorder is dynamic (a single atom oscillating between two si tes or static) or static (with a random distribution of atoms filling half of the available sites)  Figure 7 5. Structure of the high temperatur e phase of LaBGeO 5 In our calculations, the HT phase is constructed by placing atoms at the experimental positions given in Table 7 4 and generating a unit cell using the symmetry operators of the P 3 1 21 space group. Results are given as differences from neutron powder data measurements performed at 535C  To account for the half occupancy of the B and O5 sites, we make the approximation that pairs o f neighboring sites share a single atom which is placed at the mean position between the two sites. By doing so, the B and O5 atoms are forced to lie on the twofold axis to maintain the
141 P 3 1 21 symmetry as well as the stoichiometry of the cell. A schematic o f the crystal structure for the HT phase of LaBGeO 5 is shown in Fig. 7 5. In the figure, the B tetrahedra are shown in blue, Ge tetrahedra in green, La atoms in gray and O atoms in red with indices of all atoms labeled. The equilibrium lattice parameters f rom our DFT calculations using the PBEsol+ U functional are a = 6.9939 and c = 6.9508 ; LDA+ U gives a = 6.9365 and c = 6.8993 At 535C the experimental lattice parameters are a = 6.9926 and c = 6.9315  Table 7 4 Comparison of atom positions for the high temperature phase of LaBGeO 5 Experiment (535C) This work (PBEsol+ U ) Atom Position and multiplicity x/a y/b z/c x y z La 3(a) 0.4076 0 1/3 0.0067 0.0000 0.0000 B 6(c) 0.8921 0.9988 0.3226 0.0016 0.0012 0.0107 Ge 3(b) 0.4198 0 5/6 0.0005 0.0000 0.0000 O1 6(c) 0.1516 0.3418 0.0157 0.0045 0.0015 0.0004 O3 6(c) 0.1444 0.6104 0.3392 0.0034 0.0014 0.0078 O5 6(c) 0.9955 0.0367 0.8031 0.0246 0.0366 0.0303 Note: theoretical results are given as differences from experiment The fractional atom positions for the relaxed structure are given in Table 7 4 along with data from neutron powder experiments at 535C  The average deviation in our atom positions from the neutron diffraction data is < 0.1 , with the largest deviation of 0.31 for the O5 atom, which is expected as o ur approximation of the HT phase places the O5 atom between the two half occupied sites observed in experiment. In the HT phase, both the B and Ge atoms retain fourfold coordination with the surrounding oxygen, whereas the coordination of the La polyhedra increases by
142 one, see Fig. 7 6. Interatomic distances for the BO 4 and GeO 4 tetrahedra and LaO 10 polyhedra in the HT phase are given in Table 7 5 Calculated values are given as differences from neutron powder data measurements performed at 535C  Overall DFT gives good predictions for the oxygen bond distances in these polyhedra except for the La O5 bond, which is overestimated by ~0.2 Figure 7 6 (a) Structure of the LaO 9 polyhedron from the low temperature phase of LaBGeO 5 and (b) the LaO 10 polyhedron f rom the high temperature phase. 7.3.3 Lattice Dynamics of the Ferroelectric Phase Transition Calculation of the phonon dispersion of the HT phase of LaBGeO 5 reveals an unstable mode, see the dashed curve in Fig 7 7. Note ne gative numbers correspond to imaginary phonon M K near the L po int. Displacing the atoms along the eigenvectors of the unstable mode at each of the high well structures in the total energy as shown in Fig. 7 8. The double well structure for the L point modulation is omitted from the figure as the minima were much shallower than those obtained for the M K point modulations. Of these three double well structures, point eigendisplacements produce the deepest energy minima and the modulated
143 Table 7 5 Co mparison of B O, Ge O and La O interatomic distances (in ) for the high temperature phase of LaBGeO 5 Experiment (20C) This work (PBEsol+ U ) La O () O1 (x2): 2.680 2.598 O1 (x2): 0.00 2 0.031 2.771 O5 (x2): 2.747 (x2): 0.002 O5 (x2): 0.199 O3 (x2): 2.408 O3 (x2): 0.004 B O () O1: 1.516 O5: 1.451 O1: 0.014 O5: 0.007 1.516 1.451 0.014 0.007 Ge O () O1 (x2): 1.760 O1 (x2): 0.022 O3 (x2): 1.697 O3 (x2): 0.031 Note: t heoretical distances are reported as differences from experiment.
144 structures have the same symmetry as observed in the LT phase. This confirms the LT phase observed in experiments is indeed the lowest energy structure of LaBGe O 5 that can be formed from breaking the symmetry of the HT phase. Figure 7 7. Phonon dispersion and DOS for the high temperature phase of LaBGeO 5 calculated u sing DFT. is a rigid rotation of the BO 4 tetrahedra about the O1 O2) edges, which results in a z displacement of 0.16 for the O5 atoms and a much smaller z displacement of ~0.06 for the B atoms. The uncompensated z displacements of these atoms leads to a chan ge in the polarization of the cell along the c axis. This is in agreement with previous calculations based on a Valence Force Field model  that identified this mode as the main driver for the ferroelectric phase transition in LaBGeO 5 In contrast, the unstable modes at the M, K and L points are antiferrodistort ive, not producing any
145 change in polarization. In addition to the rigid BO 4 point eigendisplacements feature a twisting of the GeO 4 tetrahedra within the xy plane as well as small ~0.01 z displacements of the La and Ge atoms in the direc tion opposite to that of the B and O5 atoms. Figure 7 8. Energy as a function of the relative displacement of atoms from the high temperature phase of LaBGeO 5 and K points. By performing full relaxation calculations on the structures obtained at the minima of the double well in Fig. 7 8, we obtain two symmetry equivalent structures: P 3 1 and ( P 3 1 ) belonging to the LT phase and spatially inverted in the c direction. A CI NEB calculation between the two str uctures gives the MEP as shown in Fig. 7 9 (a) The saddle point along the MEP corresponds to a structure belonging to the high symmetry space group P 3 1 21 but with lattice parameters fixed to those of the LT phase. The height of the barrier is 0.62 eV. Thi s is a much larger energy than that of the minima of
146 the double well shown in Fig. 7 8, which indicates more substantial atomic motions in the transformation between the P 3 1 and P 3 1 21 point modulation alone. Table 7 6 giv es the atomic displacements along the MEP between the P 3 1 and P 3 1 21 point modulation from the structure at the left hand minima of the double well to the unmodulated structure (HT phase of LaBGeO 5 ). Figure 7 9. (a) MEP between the low temperature (ferroelectric) and high temperature (paraelectric) phases of LaBGeO 5 and (b) magnitude of polarization along the c axis for different points along the path shown in (a) 7.3.4 Spontaneous Polarization According to the modern theory of polarization [119 123] the spontaneous polarization P s is defined as the one ha lf the difference in polarization between two enantiomorphous polar structure to the other. The P 3 1 and ( P 3 1 ) structures calculated in the previous section correspond to such a pair of enantiomorphous structures for the LT phase of LaBGeO 5 The calculated MEP for the transformation from the P 3 1 structure to ( P 3 1 ) structure i s
147 shown in Fig. 7 9 (a). T he total polarization along the c axis calculated along the transformation path is shown in Fig. 7 9(b). As only changes i n polarization are measurable, values in Fig. 7 9(b) are given relative to the polarization at the transition state (paraelectric P 3 1 21 structure). Changes to the polarization perpendicular to the c axis are negligible. The spontaneous polarization P s is therefore directed along the c axis and is given by the difference in polarization between either of the two endpoint structures : P 3 1 and ( P 3 1 ) a nd the transition state structure P 3 1 21. For calculations using the PBEsol+ U functional, P s = 4.6 C/cm 2 whereas, without the Hubbard U correction, P s = 5.2 C/cm 2 Both calculated values lie within the experimental range of 2.7 12 C/cm 2 [98 100] 7.4 Discussion The structural electronic and vibrational properties for both the LT and HT phases of LaBGeO 5 have been calculated from first principles. Phonon dispersion calculations on the relaxed HT phase revealed a zone centered unstable mode with eigendisplacements corresponding almost exclusively to a rotation of BO 4 tetrahedra about the O1 O2) edges, see the illustration of the rotation in Fig. 7 10. Valence Force Field model calculations  had previously identified this same mode, called a rigid unit mode or RUM as exhibiting strong softening when the LT structure was chang ed towards the HT phase. Our CI NEB calculations of the MEP between the two enantiomorphous LT structures showed that the complete phase transition dynamics feature a more pronounced BO 4 rotation as well as small displacements of the La and Ge atoms, see Table 7 6 The difference in polarization between the initial and saddle points along the MEP in Fig. 7 9 is 4.6 C/cm 2 w hereas when considering only the displacements due
148 point modulation the difference in polarization between the structures at the minimum and maximum of the double well in Fig. 7 8 is 1.9 C/cm 2 The large difference in these two values can be attributed to the more pronounced rotation of the BO 4 tetrahedra that occurs along the MEP between the LT and HT phases for which the z displacements of the O5 and B atoms are more than double those o point modulation alone. Table 7 6 point modulation of the high temperature phase and those between the initial and saddle points along the MEP The spontaneous polarization estimated from the Valence Force Field model used in Ref.  is 13.8 C/cm 2 whereas our DFT calculations give a value less than half tha t. A likely source of this disagreement is an overestimation of the electronic contribution to the total macroscopic polarization, which is treated classically in the Valence Force Field model. However, according to the modern theory of polarization [119 121] the electronic polarization is a feature of the global phase of the Kohn Sham point modulation Atom x () y () z () x () y () z () La 0.0000 0.0008 0.0085 0.0780 0.0124 0.0588 B 0.0001 0.0387 0.0576 0.0878 0.0400 0.1354 Ge 0.0001 0.0142 0.0132 0.0271 0.0447 0.0639 O1 0.0101 0.0255 0.0306 0.0391 0.1144 0.1319 0.0101 0.0255 0.0307 0.0248 0.0312 0.0234 O3 0.0366 0.0099 0.0028 0.1296 0.0075 0.0728 0.0368 0.0101 0.0028 0.0319 0.0044 0.0179 O5 0.0003 0.1703 0.1598 0.1224 0.3760 0.3612
149 orbitals or, more generally, the phase of the many body elec tronic wavefunction and is, therefore, strictly a quantum mechanical phenomenon. Figure 7 10. Illustration of the rigid rotation of the BO 4 tetrahedra occurring during the ferroelectric phase transition in LaBGeO 5 7.5 Summary First principles calculations on the LT and HT phases of crystalline LaBGeO 5 give fully relaxed structures in excellent agreement with atom positions determined via neutron diffraction experiments  Phonon dispersion calculations coupled with CI NEB calculations helped to elucidate the atomic displacements leading to the ferroelectric phase transition in LaBGeO 5 The main feature of the phase transition is a rigid rotation of the BO 4 tetrahedra, which has also been recently confirmed by Valence Force Field model calculations  Our DFT calculation of the spontaneous
150 polarization in LaBGeO 5 using the PBEsol+ U functional is 4.6 C/cm 2 which is in good agreement with the e xperimental range of 2.7 12 C/cm 2 [98 100]
151 CHAPTER 8 CONCLUSIONS In Chapter 3, the intrinsic defect and self diffusion behavior in ordered phases of V 2 C was investigated using DFT calculations. Calculated activation energies for vacan cy assisted vanadium self diffusion and interstitial carbon self diffusion showed good agreement with experimental values for self diffusion of Nb in NbC x  and C in V 6 C 5  However, diffusion prefactors were underestimated in both cases. Interstitial car bon atoms were found to strongly influence vanadium migration barriers between basal planes. Migration p aths closest to the carbon interstitials were found to exhibit the highest energy barriers, while those farthest away exhibited the smallest. Moreover, differences in carbon ordering among the three V 2 C structures led to small differences in activation energies (~10 kJ/mol ) for vanadium self diffusion The work in Chapter 3 confirmed results of diffusion experiments in transition metal carbides, namely th at carbon self diffusion is several orders of magnitude faster than the metal atom self diffusion. This can be attributed to the abundance of structural vacancies on the octahedral interstitial sublattice in V 2 C, which facilitate carbon self diffusion. Pre dicted self diffusion coefficients, especially for carbon, are largely unaffected by the ordering of carbon atoms on the octahedral sublattice. Therefore, it seems reasonable to assume self diffusion in disordered V 2 C should exhibit similar activation ener gies and diffusion prefactors. Failure of DFT calculations to reproduce experimental diffusion prefactors for the case of vanadium self diffusion may be due to a lack of anharmonicity or electronic entropy in the free energy calculations. Future work inclu ding at least quasi harmonic and electronic entropy contributions to the free energy may help to determine if this is indeed the case. For carbon diffusion prefactors, the
152 underestimation is so severe that it is unlikely that more accurate calculations wil l improve the agreement much. E xperiments on carbon self diffusion in single crystals of V 2 C would help to confirm the diffusion prefactors measured in V 6 C 5  and NbC 0.868  If the results are confirmed, then other mechanisms for carbon self diffusion must b e proposed to reco ncile the large disparity observed in the diffusion prefactors predicted by DFT calculations. In Chapter, 4, the behavior of Fe, Ni, Ce, Nd and U impurities in V 2 C was investigated using DFT calculations It was found that all impurity types prefer the van adium substitutional site over the much smaller octahedral site The energy required to incorporate the impurity atom into V 2 C was shown to increase with the size of the impurity. Specifically, the impurity volumes and defect energies follow the order Fe < Ni < U < Nd < Ce Both Fe and Ni are smaller in size than the original vanadium atom while Ce, Nd and U are larger. The l arge impurities, Ce, Nd and U bind strongly to neighboring vanadium vacancies, whereas impurity vacancy binding is much weaker for the small impurities, Fe and Ni. Because of their large size, Ce, Nd and U undergo significant displacements from the ideal substitutional site upon relaxation to form defect clusters with neighboring vanadium vacancies This makes the prospect of a stand ard vacancy assisted diffusion mechanism unlikely for these impurities However, Fe and Ni remain close to the ideal substitutional site which permitted calculation of migration barriers. Compared to migration barriers for vanadium self diffusion, those o f Ni were found to be lower on average, while those of Fe were higher. The work done in Chapter 5 took the first steps to address the diffusion behavior of dilute metal impurities in V 2 C. However, the formation of strongly bound defect
153 clusters for the cas e of Ce, Nd and U seems to point to a diffusion mechanism involving multiple vanadium or carbon vacancies. Future work is necessary to determine the correct lattice diffusion mechanism for such large impurities. In the case of the Fe and Ni, migration barr iers were calculated for direct exchange of the impurity atom with a neighboring vanadium vacancy. However, calculation of the complete diffusion coefficients requires more work, and should start by considering the 8 frequency model of Ghate for diffusion of dilute impurities in an hcp lattice  In addition to lattice diffusion, future studies should also consider other me chanisms including grain boundary diffusion and diffusion along dislocations, both of which may provide fast diffusion paths for impurities. In Chapter 5 the behavior of Cu Zr MGs under shock compression was investigated using MD simulations. An elastic p lastic response was found for all compositions with the shock wave consisting of a narrow elastic precursor and a sluggishly rising plastic wave. In contrast to crystalline metals, the HEL has a strong dependence on the final shock pressure, which indicates a pressure dependent yield criterion for the simulated MGs. This agrees with shock experiments on multi component MGs that also observe d an increasing P HEL with increasing shock intensit y [57,60,200] The calculated U p U s Hugoniot (Fig. 5 7) can be broken into three regions: Region I ( U p < 0.25 km/s) nominally ela stic shock wave, r egion II (0.25 < U p < 0.75 km/s) split elastic and plastic shock waves and r egion III ( U p > 0.75 km/s) single elastic plastic shock wave. The plastic branch of the U p U s Hugoniot has approximately bilinear behavior with a gradual change in slope occurri ng near the transition from r egion II to III. The yield strength in the aftershock flow of the simulated Cu x Zr 100 x MGs was found to
154 increase with increasing Cu content over the experim ental glass forming range of 30% to 70%. Such behavior is consistent with previous MD simulations and quasistatic loading experiments. One of the major prediction s of a pressure dependent HEL needs to be verified by additional shock compression experiments as most of the available experimental data is clustered around low impact velocities. It would be particularly helpful if such experiments could be performed with sub ps time resolution of the rear surface motion, so that the narrow elastic precursor may be resolved with higher fidelity In region II, both the elastic and plastic wave s are unsteady as indicated by the slight, but continual attenuation of the elastic precursor as it propagates through the MG sample. Therefore, the shock velocities in this r egion are dependent on the length of the sample that was simulated Further simulations on longer samples might determine whether shock waves generated in region II can reach a steady state, i.e., without further attenuation of the elastic precursor It is known from MD simulations that the structure and mechanical properties of Cu x Zr 100 x MGs have a strong dependence on the cooling rate used to quench the melt [219,269 271] or, more generally, on the configurational potential energy of the resulting glassy structure  In addition quasistatic loading experiments have found that the mechanical properties of Cu x Zr 100 x MGs undergo significant changes when small amounts of Al or Ag are added [273 276] It would be interesting to see how changes in configurational potential energy or small element additions affect the shock response of the resulta nt MGs. Such studies would be facilitated by the availability of interatomic potentials for both the Cu Zr Al  and Cu Zr Ag  systems that were both fit to amorphous structures.
155 In Chapter 6, the focus of the analysis was turned to the micro structural response of Cu x Zr 100 x MGs during shock compression It was found that plastic deformation occurs via nucleation of localized regions of higher than average shear strain, termed shear transformation zones (S TZs ) With increasing shock intensity, the rate of STZ nucleation behind the shock front increase d with a gradual progression towards homogeneous plastic flow. Moreover, the amount of shear strain at a given shock pressure wa s also shown to decrease with increasing Cu content which is consistent with the increase in flow stress with Cu content observed in Chapter 5 T he transition to shock induced melting was also found to be dependent on the MG composition with higher Cu content leading to higher shock pressures for melting. An alysis of the local structure showed that regions of h igh icosahedral order (<0,0,12,0> Voronoi polyhedra) were most resist ant to shear deformation As Cu x Zr 100 x MGs with higher Cu composition s exhibit more icosahedral order, it is clear why we observe the similar trends in the composition dependence of the flow strength, shear strain and melting transition discussed above T he shock pressure at the melting transition found in Chapter 6 was estimated using the diffusivity. From this analysis alone, it was not clear how the thermodynamic state of the material with in the aftershock flow compared to the equilibrium liquid at the same pressure. By comparing the temperature and pressure of the material in the aftershock flow to the equilibrium melting line T m ( P ), one could dete rmine whether the material is superheated, i.e., a solid above the equilibrium melting line, or sup ercooled i.e., a liquid below the equilibrium melting line. As was shown in previous MD simulations of single crystal Al, recrystallization can occur behind the shock front if the
156 melt is sufficiently undercooled  Further work on shock induced melting of Cu x Zr 100 x MGs should calculate the theoretical equilibrium melting line for different compositions. A way of differentiating solid amorphous structures from liquid ones also needs to be developed to complement the diffusivity analysis used in Chapter 6. One possible method is an in plane RDF, which was used succes sfully in Ref.  to differentiate crystalline solid and liq uid atomic configurations. It has be en shown that MGs exhibit ordering that extends a few nanometers beyond the first nearest neighbor shell [278 280] In Cu x Zr 100 x MGs, this so ca lled medium range order develops from connections among Cu centered clusters typically clusters of icosahedra [91,245,259,281,282] MD simulations have found th at [283 287] Future work on shock loading of Cu x Zr 100 x MGs should focus on the evolution of connections among Cu centered clusters as a function of shock intensity. One method to quantitatively assess the connectivity of d ifferent clusters is through the strength of the spatial correlation for different Voronoi polyhedra types [222,282] In Chapter 7, the dynamics of the ferroelectric phase transition was investigated in the stillwellite compound LaBGeO 5 From DFT c alculations of the phonon dispersion for the high temperature paraelectric phase it was fou nd that an unstable mode exists Displacing the atoms along the eigenvect ors of this mode led to two equivalent structures with oppositely oriented polarization vectors that had the symmetry of the low temperature ferroelectric phase. The main feature of the displace ments is a rigid rotation of BO 4 tetrahedra, which is consiste nt with previous experimental and theoretical results [117,118] The spontaneous polarization was
157 determined by calculating the change in the total p olarization along a continuous path between the two enantiomorphous polar structures. The calculated value of 4.6 C/cm 2 lie s within the experimental range of 2.7 12 C/cm 2 [98 100] In this work, we used atomic scale simulation methods to explore the thermodynamics and kinetics of defect formation and diffusion, phase transitions and high strain rate plastic deformation processes in ceramics and metals. Both static and dynamic simulation methods were used to explore different aspects of the PES. Energy and force based energy minimization methods were used to determine local energy minima to predict various defect structures and crystalline phases for V 2 C and LaBGeO 5 Constrained search methods, like the CI NEB method, were used to determine MEPs and, hence, transition states for diffusive hops in V 2 C. Force constant matrix calculations of the transition states in V 2 C and LaBGeO 5 were used to determine eigendisplacements and phonon frequencies. The process of glass formation in the Cu Zr system was simulated using MD through rapid cooling of the melt to a quasi stable local energy minimum corresponding to the amorphous state. Non equilibrium MD landscape in response to high strain rate shock loading. Small clusters of atoms were shown to exhibit thermally activated shear deformation events under suff iciently high stresses. Under action of a strong shock, the MG system was shown to undergo melting as it collectively escapes the local energy basin of the glassy state.
158 LIST OF REFERENCES  P. G. Debenedetti and F. H. Stillinger, Nature 410 259 (2001).  C. A. Schuh, T. C. Hufnagel, and U. Ramamurty, Acta Mater. 55 4067 (2007).  G. L. Hofman, A. G. Hins, D. L. Porter, L. Leibowitz, and E. L. Wood, Int. Conf. Reliab. Fuels Liq. Met. Rea ct. (1986).  L. A. Lawrence, Nucl. Technol. 64 139 (1984).  C. Gueneau, J. P. Piron, J. C. Dumas, V. Bouineau, F. C. Iglesias, and B. J. Lewis, Fuel Cladding Chemical Interaction (Nuclear Energy Agency of the OECD (NEA), 2015).  D. C. Fee and C. E. Johnson, J. Nucl. Mater. 96 80 (1981).  R. J. M. Konings, T. R. Allen, R. E. Stoller, and S. Yamanaka, Comprehensive Nuclear Materials (Elsevier Science, Burlington, 2011).  D. D. Keiser and J. I. Cole, TMS Lett. 2 79 (2005).  D. D. Keiser an d M. A. Dayananda, J. Nucl. Mater. 200 229 (1993).  R. G. Pahl, C. E. Lahm, and S. L. Hayes, J. Nucl. Mater. 204 141 (1993).  A. B. Cohen, H. Tsai, and L. A. Neimark, J. Nucl. Mater. 204 244 (1993).  H. J. Ryu, B. O. Lee, S. J. Oh, J. H. Kim and C. B. Lee, J. Nucl. Mater. 392 206 (2009).  S. W. Yang, H. J. Ryu, J. H. Kim, B. O. Lee, and C. B. Lee, J. Nucl. Mater. 401 98 (2011).  D. D. Keiser and J. I. Cole, An Evaluation of Potential Liner Materials for Eliminating FCCI in Irradiat ed Metallic Nuclear Fuel Elements (Idaho National Laboratory (INL), 2007).  W. Y. Lo and Y. Yang, J. Nucl. Mater. 451 137 (2014).  S. Sarian, J. Appl. Phys. 39 3305 (1968).  F. J. J. Van Loo, W. Wakelkamp, G. F. Bastin, and R. Metselaar, Solid State Ionics 32 824 (1989).  S. Sarian and J. M. Criscione, J. Appl. Phys. 38 1794 (1967).  B. B. Yu and R. F. Davis, Phys. Status Solidi 51 261 (1979).
159  B. B. Yu and R. F Davis, J. Phys. Chem. Solids 40 997 (1979).  A. V Shovensin, G. V Shcherbedinskii, and A. N. Minkevich, Sov. Powder Metall. Met. Ceram. 5 880 (1966).  S. Sarian, J. Appl. Phys. 39 5036 (1968).  W. Brizes, L. Cadoff, and J. Tobin, J. Nucl. Mater. 20 57 (1966).  S. Sarian, J. Phys. Chem. Solids 33 1637 (1972).  C. P. Buhsmer and P. H. Crayton, J. Mater. Sci. 6 981 (1971).  G. L. DePoorter and T. Wallace, Adv. High Temp. Chem. 4 107 (1970).  B. B. Yu and R. F. Davis, J. Phy s. Chem. Solids 42 83 (1981).  S. Sarian, J. Appl. Phys. 40 3515 (1969).  M. Mantina, Mater. Sci. Eng. PhD 252 (2008).  M. Mantina, Y. Wang, R. Arroyave, L. Q. Chen, Z. K. Liu, and C. Wolverton, Phys. Rev. Lett. 100 1 (2008).  M. Mantina, Y. Wang, R. Arroyave, S. L. Shang, L. Q. Chen, and Z. K. Liu, J. Phys. Condens. Matter 24 305402 (2012).  M. Mantina, Y. Wang, L. Q. Chen, Z. K. Liu, and C. Wolverton, Acta Mater. 57 4102 (2009).  M. Mantina, L. Q. Chen, and Z. K. Liu Defect Diffus. Forum 294 1 (2009).  S. Huang, D. L. Worthington, M. Asta, V. Ozolins, G. Ghosh, and P. K. Liaw, Acta Mater. 58 1982 (2010).  N. Sandberg and R. Holmestad, Phys. Rev. B Condens. Matter Mater. Phys. 73 1 (2006).  S. Ganesha n, L. G. Hector, and Z. K. Liu, Acta Mater. 59 3214 (2011).  S. Ganeshan, L. G. Hector, and Z. K. Liu, Comput. Mater. Sci. 50 301 (2010).  C. Z. Hargather, S. L. Shang, Z. K. Liu, and Y. Du, Comput. Mater. Sci. 86 (2014).  M. Chen, Annu. Rev. Mater. Res. 38 445 (2008).  M. D. Demetriou, M. E. Launey, G. Garrett, J. P. Schramm, D. C. Hofmann, W. L. Johnson, and R. O. Ritchie, Nat. Mater. 10 123 (2011).
160  A. L. Greer, Mater. Today 12 14 (2009).  A. L. Greer and E Ma, MRS Bull. 32 611 (2007).  W. L. Johnson, JOM 54 40 (2002).  J. F. Lffler, Intermetallics 11 529 (2003).  M. Miller and P. Liaw, Bulk Metallic Glasses (2008).  C. Suryanarayana and A. Inoue, Bulk Metallic Glasses (CRC Press, 2010).  M. M. Trexler and N. N. Thadhani, Prog. Mater. Sci. 55 759 (2010).  W. H. Wang, C. Dong, and C. H. Shek, Mater. Sci. Eng. R Reports 44 45 (2004).  A. R. Yavari, J. J. Lewandowski, and J. Eckert, MRS Bull. 32 635 (2007).  T. C. Hufnagel C. A. Schuh, and M. L. Falk, Acta Mater. 109 375 (2016).  T. Mukai, T. G. Nieh, Y. Kawamura, A. Inoue, and K. Higashi, Intermetallics 10 1071 (2002).  T. C. Hufnagel, T. Jiao, Y. Li, L. Q. Xing, and K. T. Ramesh, J. Mater. Res. 17 1441 (2002).  R. F. Trunin, Shock Compression of Condensed Materials (Cambridge University Press, Cambridge, 1998).  L. Davison and R. A. Graham, Phys. Rep. 55 255 (1979).  L. C. Chhabildas, L. Davison, and Y. Horie, High Pressure Shock Compression of Solids VIII: The Science and Technology of High Velocity Impact (Springer Berlin Heidelberg, 2004).  S. J. Turneaure, J. M. Winey, and Y. M. Gupta, Appl. Phys. Lett. 84 1692 (2004).  S. J. Turneaure, J. M. Winey, and Y. M. Gupta, J. Appl. Phys. 10 0 63522 (2006).  T. Mashimo, H. Togo, Y. Zhang, Y. Uemura, T. Kinoshita, M. Kodama, and Y. Kawamura, Appl. Phys. Lett. 89 31 (2006).  H. Togo, Y. Zhang, Y. Kawamura, and T. Mashimo, Mater. Sci. Eng. A 448 451 264 (2007).  M. Martin, T. Sekine, T. Kobayashi, L. Kecskes, and N. N. Thadhani, Metall. Mater. Trans. A Phys. Metall. Mater. Sci. 38 A 2689 (2007).  F. Yuan, V. Prakash, and J. J. Lewandowski, Mech. Mater. 41 886 (2009).
161  S. A. Atroshenko, N. F. Morozov, W. Zheng, Y. J. Huang, Y. V. Sudenkov, N. S. Naumova, and J. Shen, J. Alloys Compd. 505 501 (2010).  F. Xi, Y. Yu, C. Dai, Y. Zhang, and L. Cai, J. Appl. Phys. 108 1 (2010).  T. Jaglinski, S. J. Turneaure, and Y. M. Gupta, J. Appl. Phys. 112 1 (201 2).  B. Luo, G. Wang, F. Tan, J. Zhao, C. Liu, and C. Sun, AIP Adv. 5 (2015).  Y. Y. Yu, F. Xi, C. Da Dai, L. C. Cai, Y. Tan, X. M. Li, Q. Wu, and H. Tan, Chinese Phys. B 24 66201 (2015).  Y. Hironaka, A. Yazaki, F. Saito, K. G. Nakamura, K. Kondo, H. Takenaka, and M. Yoshida, Appl. Phys. Lett. 77 1967 (2000).  S. N. Luo, D. C. Swift, T. E. Tierney, D. L. Paisley, G. A. Kyrala, R. P. Johnson, A. A. Hauer, O. Tschauner, and P. D. Asimow, High Press. Res. 24 409 (2004).  K. Rosolankova J. S. Wark, E. M. Bringa, and J. Hawreliak, J. Phys. Condens. Matter 18 6749 (2006).  M. Righini, R. Quidant, S. R. Emory, J. A. Dionne, H. A. Atwater, Y. Kadoya, H. F. Hofmann, Q. Wang, F. Capasso, V. M. Shalaev, A. Polman, A. Archambault, T. V Tep erik, F. Marquier, J. J. Greffet, P. Yeh, D. Milathianaki, S. Boutet, G. J. Williams, A. Higginbotham, D. Ratner, A. E. Gleason, M. Messerschmidt, M. M. Seibert, D. C. Swift, P. Hering, J. Robinson, W. E. White, and J. S. Wark, Science (80 ). 342 220 (2 013).  M. G. Gorman, R. Briggs, E. E. McBride, A. Higginbotham, B. Arnold, J. H. Eggert, D. E. Fratanduono, E. Galtier, A. E. Lazicki, H. J. Lee, H. P. Liermann, B. Nagler, A. Rothkirch, R. F. Smith, D. C. Swift, G. W. Collins, J. S. Wark, and M. I. Mc Mahon, Phys. Rev. Lett. 115 1 (2015).  G. Mogni, A. Higginbotham, K. Gal Nagy, N. Park, and J. S. Wark, Phys. Rev. B 89 64104 (2014).  T. R. Mattsson, J. M. D. Lane, K. R. Cochrane, M. P. Desjarlais, A. P. Thompson, F. Pierce, and G. S. Grest, P hys. Rev. B 81 54103 (2010).  J. Budzien, A. P. Thompson, and S. V. Zybin, J. Phys. Chem. B 113 13142 (2009).  H. N. Jarmakani, E. M. Bringa, P. Erhart, B. A. Remington, Y. M. Wang, N. Q. Vo, and M. A. Meyers, Acta Mater. 56 5584 (2008).  M. J. Cawkwell, T. D. Sewell, L. Zheng, and D. L. Thompson, Phys. Rev. B Condens. Matter Mater. Phys. 78 1 (2008).
162  S. Zhao, T. C. Germann, and A. Strachan, Phys. Rev. B Condens. Matter Mater. Phys. 76 1 (2007).  T. Germann, B. Holian, P. Lomda hl, and R. Ravelo, Phys. Rev. Lett. 84 5351 (2000).  83 1175 (1999).  D. H. Robertson, D. W. Brenner, and C. T. White, Phys. Rev. Lett. 67 3132 (1991).  E. M. Bringa, A. Caro, Y. Wang, M. Victoria, J. M. McNaney, B. A. Remington, R. F. Smith, B. R. Torralva, and H. Van Swygenhoven, Science 309 1838 (2005).  K. Kadau, T. C. Germann, P. S. Lomdahl, and B. L. Holian, Phys. Rev. B Condens. Matter Mater. Phys. 72 1 (2005).  E. M. Bringa, K. Rosolankova, R. E. Rudd, B. A. Remington, J. S. Wark, M. Duchaineau, D. H. Kalantar, J. Hawreliak, and J. Belak, Nat. Mater. 5 805 (2006).  H. N. Jarmakani, E. M. Bringa, P. Erhart, B. A. Remington, Y. M. Wang, N. Q. Vo, and M. A. Meyers, Acta Mater. 56 5584 (2008).  V. V. Zhakhovsky, M. M. Budzevich, N. A. Inogamov, I. I. Oleynik, and C. T. White, Phys. Rev. Le tt. 107 1 (2011).  M. M. Budzevich, V. V. Zhakhovsky, C. T. White, and I. I. Oleynik, Phys. Rev. Lett. 109 (2012).  W. Ma, W. Zhu, and Y. Hou, J. Appl. Phys. 114 (2013).  R. Ravelo, T. C. Germann, O. Guerrero, Q. An, and B. L. Holian, Phys. Rev. B Condens. Matter Mater. Phys. 88 (2013).  D. Tramontina, P. Erhart, T. Germann, J. Hawreliak, A. Higginbotham, N. Park, R. Ravelo, A. Stukowski, M. Suggit, Y. Tang, J. Wark, and E. Bringa, High Energy Density Phys. (2014).  K. Chen, Y. Yu, Z. Zhang, and S. Q. Shi, Scr. Mater. 120 62 (2016).  W. R. Jian, X. H. Yao, L. Wang, X. C. Tang, and S. N. Luo, J. Appl. Phys. 118 15901 (2015).  P. Wen, G. Tao, C. Pang, S. Yuan, and Q. Wang, Comput. Mater. Sci. 124 304 (2016).  B. Arman, Condens. Matter Mater. Phys. 81 1 (2010).
163  E. M. Bringa, Science (80 ). 309 1838 (2005).  V. N. Sigaev, E. V. Lopatina, P. D. Sarkisov, S. Y. Stefanovich, and V. I. Molev, Mater. Sci. Eng. B 48 254 (1997).  V. N. Sigaev, S. Y. Stefanovich, P. D. Sarkisov, and E. V Lopatina, Glas. Phys. Chem. 20 392 (1994).  V. N. Sigaev, S. Y. Stefanovich, P. D. Sarkisov, and E. V. Lopatina, Mater. Sci. Eng. B 32 17 (1995).  S. Y. Stefanovich, B. V Mill, and A. V Butashin, Kristallografiya 37 965 (1992).  A. Onodera, B. A. Strukov, and A. Belov, J. Phys. Soc. Japan 62 4311 (1993).  E. V. Milov and B. A. Strukov, Phys. Solid State 43 513 (2001).  E. V. Milov, B. A. Strukov, and V. N. Milov, Ferroelectrics 269 15 (2002).  B. A. Strukov, E. V. Milov, V. N. Milov, A. P. Korobtsov, T. Tomida, K. Sato, M. Fukunaga, and Y. Uesu, Ferroelectrics 314 105 (2005). [ 103 ] J. Capmany, D. Jaque, J. Garc a Sol, and A. A. Kaminskii, Appl. Phys. Lett. 72 531 (1998). [ 104 ] J. Capmany and J. Garc a Sol, Appl. Phys. Lett. 70 2517 (1997).  Y. Takahashi, Y. Benino, V. Dimitrov, and T. Komatsu, J. Non. Cryst. Solids 260 155 (1999).  Y. Takahashi, Y. Benino, T. Fujiwara, and T. Komatsu, J. Appl. Phys. 89 5282 (2001).  V. N. Sigaev, P. D. Sarkisov, P. Pernice, A. Aronne, A. M. Datsenko, S. Y. Stefanovich, V. I. Fertikov, O. A. Pozhogin, and D. A. Zakharkin, J. Eur. Ceram. Soc. 24 1063 (2004). [1 08] P. Gupta, H. Jain, D. B. Williams, O. Kanert, and R. Kuechler, J. Non. Cryst. Solids 349 291 (2004).  V. N. Sigaev, S. V. Lotarev, E. V. Orlova, S. Y. Stefanovich, P. Pernice, A. Aronne, E. Fanelli, and I. Gregora, J. Non. Cryst. Solids 353 1956 (2007).  V. N. Sigaev, E. A. Alieva, S. V. Lotarev, N. M. Lepekhin, Y. S. Priseko, and A. V. Rasstanaev, Glas. Phys. Chem. 35 13 (2009).  A. Stone, M. Sakakura, Y. Shimotsuma, G. Stone, P. Gupta, K. Miura, K. Hirao, V. Dierolf, and H. Jain, Opt Express 17 23284 (2009).
164  Proc. SPIE (2015), p. 945018.  A. Stone, H. Jain, V. Dierolf, M. Sakakura, Y. Shimotsuma, K. Miura, K. Hirao, J. Lapointe, and R. Ka shyap, Sci. Rep. 5 10391 (2015).  N. Sigaev, Glas. Ceram. (English Transl. Steklo I Keramika) 73 441 (2017).  G. Okhrimchuk, A. S. Larkin, M. Y. Presnyakov, and V. N. Sigaev, Cryst. Growth Des. 17 4670 (2017).  E. L. Belokoneva, W. I. F. David, J. B. Forsyth, and K. S. Knight, J. Physics. Condens. Matter 9 3503 (1997).  S. Kamba, J. Petzelt, I. Grego ra, and Z. Zikmund, Phys. Stat. Sol. 214 (1999).  M. B. Smirnov, A. V. Menschikova, I. Kratochvilova Hruba, and Z. Zikmund, Phys. Status Solidi 241 1017 (2004).  R. D. King Smith and D. Vanderbilt, Phys. Rev. B 47 1651 (1993).  D. Vanderb ilt and R. D. King Smith, Phys. Rev. B 48 4442 (1993).  R. Resta, Rev. Mod. Phys. 66 899 (1994).  R. Resta, Model. Simul. Mater. Sci. Eng. 11 R69 (2003).  R. Resta and D. Vanderbilt, Top. Appl. Phys. 105 31 (2007).  G. H. Vineyard, J. Phys. Chem. Sol. 3 121 (1957).  M. Mantina, A First Principles Methodology for Diffusion Coefficients in Metals and Dilute Alloys (ProQuest, 2008).  P. Hohenberg and W. Kohn, Phys. Rev. 136 B864 (1964).  W. Kohn and L. J Sham, Phys. Rev. 140 A1133 (1965).  J. P. Perdew and A. Zunger, Phys. Rev. B 23 5048 (1981).  J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77 3865 (1996).  J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scu seria, L. A. Constantin, X. Zhou, and K. Burke, Phys. Rev. Lett. 100 136406 (2008).  G. Kresse and J. Furthmller, Phys. Rev. B 54 11169 (1996).  N. A. Spaldin, J. Solid State Chem. 195 2 (2012).
165  L. Verlet, Phys. Rev. 159 98 (1967). [13 4] W. C. Swope, H. C. Andersen, P. H. Berens, and K. R. Wilson, J. Chem. Phys. 76 637 (1982).  C. W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations (Prentice Hall, 1971).  S. Nos, J. Chem. Phys. 81 511 (1984).  S. Nos, Mol. Phys. 52 255 (1984).  W. G. Hoover, Phys. Rev. A 31 1695 (1985).  M. Parrinello and A. Rahman, J. Appl. Phys. 52 7182 (1981).  S. Plimpton, J. Comput. Phys. 117 1 (1995).  J. E. Jones, Proc. R. Soc. London. Ser. A 106 463 LP (1924).  A. C. T. van Duin, S. Dasgupta, F. Lorant, and W. A. Goddard, J. Phys. Chem. A 105 9396 (2001).  M. S. Daw and M. I. Baskes, Phys. Rev. B 29 6443 (1984).  M. S. Daw, S. M. Foiles, and M. I. Baskes, Mater. Sci. Reports 9 2 51 (1993).  Y. Mishin, in Handb. Mater. Model. SE 23 edited by S. Yip (Springer Netherlands, 2005), pp. 459 478.  L. E. Toth, Transition Metal Carbides and Nitrides (Academic Press, New York, 1971).  S. T. Oyama, The Chemistry of Transiti on Metal Carbides and Nitrides (Springer Science & Business Media, 1996).  H. O. H. Pierson, Handb. Refract. Carbides Nitrides Prop. Charact. Process. Apps. 81 (1996).  D. D. Keiser and J. Cole, Mater. Sci. 597 (2007).  O. N. Carlson, A. H. Ghaneya, and J. F. Smith, Bull. Alloy Phase Diagrams 6 115 (1985).  V. N. Lipatnikov, A. I. Gusev, P. Ettmayer, and W. Lengauer, J. Phys. Condens. Matter 11 163 (1999).  V. N. Lipatnikov, Russ. Chem. Rev. 74 697 (2005).
166  C. P. Buhsmer and P. H. Crayton, J. Mater. Sci. 6 981 (1971).  B. Ozturk, V. L. Fearing, J. A. Ruth, and G. Simkovich, Metall. Trans. A 13 1871 (1982).  X. Chong, Y. Jiang, R. Zhou, and J. Feng, RSC Adv. 4 44959 (2014).  L. Wu, T. Yao, Y. Wang, J. Zhang, F. Xiao, and B. Liao, J. Alloys Compd. 548 60 (2013).  H. Liu, J. Zhu, Y. Liu, and Z. Lai, Mater. Lett. 62 3084 (2008).  H. Matzke, in Defect Diffus. Forum (Trans Tech Publ, 1992), pp. 111 130.  G. Kresse and J. Furthmller, Comput. Mat. S ci. 6 15 (1996).  G. Kresse and D. Joubert, Phys. Rev. B 59 1758 (1999).  P. E. Blchl, Phys. Rev. B 50 17953 (1994).  P. Pulay, Chem. Phys. Lett. 73 393 (1980).  P. E. Blchl, O. Jepsen, and O. K. Andersen, Phys. Rev. B 49 16223 (1994).  A. Togo, F. Oba, and I. Tanaka, Phys. Rev. B Condens. Matter Mater. Phys. 78 1 (2008).  X. Wu, D. Vanderbilt, and D. R. Hamann, Phys. Rev. B Condens. Matter Mater. Phys. 72 1 (2005).  A. U. Ortiz, A. Boutin, A. H. Fuchs, and F X. Coudert, J. Chem. Phys. 138 174703 (2013).  A. Togo and I. Tanaka, Scr. Mater. 108 1 (2015).  G. Mills, H. Jnsson, and G. K. Schenter, Surf. Sci. 324 305 (1995).  H. Jnsson, G. Mills, and K. W. Jacobsen, in Class. Quantum Dyn. Cond ens. Phase Simulations (Citeseer, 1998), pp. 385 404.  E. Rudy and C. E. Brukl, J. Am. Ceram. Soc. 50 265 (1967).  S. Curtarolo, W. Setyawan, G. L. W. Hart, M. Jahnatek, R. V. Chepulskii, R. H. Taylor, S. Wang, J. Xue, K. Yang, O. Levy, M. J. Me hl, H. T. Stokes, D. O. Demchenko, and D. Morgan, Comput. Mater. Sci. 58 218 (2012).  G. L. W. Hart and R. W. Forcade, Phys. Rev. B 80 14120 (2009).
167  A. K. Rappe, C. J. Casewit, K. S. Colwell, W. A. Goddard, and W. M. Skiff, J. Am. Chem. Soc. 114 10024 (1992).  W. Xing, F. Meng, and R. Yu, Sci. Rep. 6 21794 (2016).  M. P. Arbuzov, V. G. Fak, and B. V Khaenko, Izv. Akad. Nauk SSSR, Neorg. Mater 12 846 (1976).  K. Yvon, H. Nowotny, and R. Kieffer, Monatshefte Fr Chemie/Chemical Mon. 98 34 (1967).  R. Hill, Proc. Phys. Soc. Sect. A 65 349 (2002).  N. W. Ashcroft and N. D. Mermin, Solid State Physics (Saunders College, 1976).  S. Lany and A. Zunger, Model. Simul. Mater. Sci. Eng. 17 84002 (2009).  C. Freysoldt, B. Grabowski, T. Hickel, J. Neugebauer, G. Kresse, A. Janotti, and C. G. Van De Walle, Rev. Mod. Phys. 86 253 (2014).  M. Hagen, Philos. Mag. A 77 447 (1998).  C. G. Van De Walle and J. Neugebauer, J. Appl. Phys. 95 3851 (2004).  M. T. Dove, Introduction to Lattice Dynamics (Cambridge University Press, 1993).  J. G. Mullen, Phys. Rev. 124 1723 (1961).  L. V. Azroff, J. Appl. Phys. 32 1658 (1961).  L. V. Azroff, J. Appl. Phys. 32 1663 (1961).  J. E. Sa al and C. Wolverton, Acta Mater. 60 5151 (2012).  S. L. Dudarev, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, Phys. Rev. B 57 1505 (1998).  W. Xie, W. Xiong, C. A. Marianetti, and D. Morgan, Phys. Rev. B 88 235128 (2013).  M. Methfessel and A. T. Paxton, Phys. Rev. B 40 3616 (1989).  W. M. Haynes, D. R. Lide, and T. J. Bruno, CRC Handbook of Chemistry and Reference Book of Chemical and Physical Data. (2015).  G. Henkelman, B. P. Uberuaga, and H. Jnsson, J. Chem. Phys. 113 9901 (2000).  D. Sheppard, R. Terrell, and G. Henkelman, J. Chem. Phys. 128 134106 (2008).
168  B. J. Demaske, A. Chernatynskiy, and S. R. Phillpot, J. Phys. Condens. Matter 29 (2017).  D. Shin and C. Wolverton, Ac ta Mater. 58 531 (2010).  C. Wolverton, Acta Mater. 55 5867 (2007).  A. Inoue and A. Takeuchi, Acta Mater. 59 2243 (2011).  Y. Q. Cheng and E. Ma, Prog. Mater. Sci. 56 379 (2011).  W. H. Wang, Prog. Mater. Sci. 57 487 (2012).  F. Xi, Y. Yu, C. Dai, Y. Zhang, and L. Cai, J. Appl. Phys. 108 1 (2010).  X. Huang, Z. Ling, H. S. Zhang, J. Ma, and L. H. Dai, J. Appl. Phys. 110 103519 (2011).  L. Lu, C. Li, W. H. Wang, M. H. Zhu, X. L. Gong, and S. N. Luo, Mater. Sci. Eng. A 651 848 (2016).  B. Cao, E. M. Bringa, and M. A. Meyers, Metall. Mater. Trans. A 38 2681 (2007).  E. M. Bringa, J. U. Cazamias, P. Erhart, J. Stlken, N. Tanushe v, B. D. Wirth, R. E. Rudd, and M. J. Caturla, J. Appl. Phys. 96 3793 (2004).  M. M. Sichani and D. E. Spearot, J. Appl. Phys. 120 45902 (2016).  M. M. Sichani and D. E. Spearot, Comput. Mater. Sci. 108 226 (2015).  W. Ma, W. Zhu, and Y. Hou, J. Appl. Phys. 114 1 (2013).  W. Ma, W. Zhu, and J. Fuqian, Appl. Phys. Lett. 97 121903 (2010).  S. Marinier and L. J. Lewis, Phys. Rev. B 92 184108 (2015).  K. G. Chen, Y. Y. Yu, Z. G. Zhang, and S. Q. Shi, Scr. Mater. 120 62 (2016 ).  X. Huang, Z. Ling, and L. H. Dai, J. Appl. Phys. 116 (2014).  S. G. Hao, C. Z. Wang, M. J. Kramer, and K. M. Ho, J. Appl. Phys. (2010).  H. L. Peng, M. Z. Li, W. H. Wang, C. Wang, and K. M. Ho, Appl. Phys. Lett. 96 30 (2010).  Z. D. Sha, Y. P. Feng, and Y. Li, Appl. Phys. Lett. 96 61903 (2010).  J. C. Lee, K. W. Park, K. H. Kim, E. Fleury, B. J. Lee, M. Wakeda, and Y. Shibutani, J. Mater. Res. 22 3087 (2007).
169  K. W. Park, J. I. Jang, M. Wakeda, Y. Shibutani, and J. C. Lee, Scr. Mater. 57 805 (2007).  Y. Q. Cheng, H. W. Sheng, and E. Ma, Phys. Rev. B Condens. Matter Mater. Phys. 78 1 (2008).  N. Mattern, A. Schps, U. Khn, J. Acker, O. Khvostikova, and J. Eckert, J. Non. Cryst. Solids 354 1054 (2008).  Y. Q. Cheng, A. J. Cao, H. W. Sheng, and E. Ma, Acta Mater. 56 5263 (2008).  N. Mattern, P. Jovari, I. Kaban, S. Gruner, A. Elsner, V. Kokotin, H. Franz, B. Beuneu, and J. Eckert, J. Alloy Compd. J. 485 163 (2009).  A. E. Lagogianni, G. Almyras, C. E. Lekka, and D. G. Papageorgiou, J. Alloys Compd. 483 658 (2009).  M. Li, C. Wang, S. Hao, M. Kramer, and K. Ho, Phys. Rev. B 80 1 (2009).  M. I. Mendelev, D. J. Sordelet, and M. J. Kramer, J. Appl. Phys. 102 (2007).  Y. Q. Cheng, E. Ma, and H. W. Sheng, Phys. Rev. Lett. 102 (2009).  L. Davison, Fundamentals of Shock Wave Propagation in Solids (Springer Berlin Heidelberg, 2008). [2 26] H. Okamoto, J. Phase Equilibria Diffus. 29 204 (2008).  B. L. Holian, Science (80 ). 280 2085 (1998).  M. A. Meyers, H. Jarmakani, E. M. Bringa, and B. A. Remington, Dislocations in Shock Compression (2009).  L. A. Davis and S. Kaves h, J. Mater. Sci. 10 453 (1975).  P. E. Donovan, Acta Metall. 37 445 (1989).  K. M. Flores and R. H. Dauskardt, Acta Mater. 49 2527 (2001).  J. J. Lewandowski and P. Lowhaphandu, Philos. Mag. A 82 3427 (2002).  P. Lowhaphandu, S. L. Montgomery, and J. J. Lewandowski, Scr. Mater. 41 19 (1999).  J. Lu and G. Ravichandran, J. Mater. Res. 18 2039 (2003).  R. Vaidyanathan, M. Dao, G. Ravichandran, and S. Suresh, Acta Mater. 49 3781 (2001).
170  A. Peker and W. L. Johnson, Ap pl. Phys. Lett. 63 2342 (1993).  O. J. Kwon, Y. C. Kim, K. B. Kim, Y. K. Lee, and E. Fleury, Met. Mater. Int. 12 207 (2006).  N. Mattern, P. Jvri, I. Kaban, S. Gruner, A. Elsner, V. Kokotin, H. Franz, B. Beuneu, and J. Eckert, J. Alloys Compd 485 163 (2009).  J. Antonowicz, A. Pietnoczka, W. Zalewski, R. Bacewicz, M. Stoica, K. Georgarakis, and A. R. Yavari, J. Alloys Compd. 509 S34 (2011).  M. I. Mendelev, D. K. Rehbein, R. T. Ott, M. J. Kramer, and D. J. Sordelet, J. Appl. Phys. 102 (2007).  M. I. Mendelev, M. J. Kramer, R. T. Ott, D. J. Sordelet, D. Yagodin, and P. Popel, Philos. Mag. 89 967 (2009).  A. S. Argon, Acta Metall. 27 47 (1979).  C. A. Schuh and A. C. Lund, Nat. Mater. 2 449 (2003).  A. C. Lund and C. A. Schuh, Acta Mater. 51 5399 (2003).  Y. Q. Cheng and E. Ma, Prog. Mater. Sci. 56 379 (2011).  K. Yoshimoto, T. S. Jain, K. Van Workum, P. F. Nealey, and J. J. de Pablo, Phys. Rev. Lett. 93 175501 (2004).  M. Tsamados, A. Tanguy, C. Goldenberg, and J. L. Barrat, Phys. Rev. E 80 26112 (2009).  S. G. Mayr, Phys. Rev. B 79 60201 (2009).  H. L. Peng, M. Z. Li, and W. H. Wang, Phys. Rev. Lett. 106 1 (2011).  F. Shimizu, S. Ogata, and J. Li, Mater. Trans. 48 2923 (2007).  C. H. Rycroft, Chaos An Interdiscip. J. Nonlinear Sci. 19 41111 (2009).  H. W. Sheng, W. K. Luo, F. M. Alamgir, J. M. Bai, and E. Ma, Nature 439 419 (2006).  A. Stukowski, Model. Simul. Mater. Sci. Eng. 18 15012 (2010). [25 4] M. Wakeda, Y. Shibutani, S. Ogata, and J. Park, Intermetallics 15 139 (2007).  C. E. Lekka, A. Ibenskas, A. R. Yavari, and G. A. Evangelakis, Appl. Phys. Lett. 91 214103 (2007).
171  Y. Q. Cheng, A. J. Cao, and E. Ma, Acta Mater. 57 3253 (2009)  C. Zhong, H. Zhang, Q. P. Cao, X. D. Wang, D. X. Zhang, U. Ramamurty, and J. Z. Jiang, Scr. Mater. 114 93 (2016).  D. Ma, A. D. Stoica, X. L. Wang, Z. P. Lu, M. Xu, and M. Kramer, Phys. Rev. B Condens. Matter Mater. Phys. 80 1 (2009). [259 ] Y. Q. Cheng, E. Ma, and H. W. Sheng, Phys. Rev. Lett. 102 1 (2009).  L. Zhang, Y. Cheng, A. Cao, J. Xu, and E. Ma, Acta Mater. 57 1154 (2009).  S. Trady, M. Mazroui, A. Hasnaoui, and K. Saadouni, J. Non. Cryst. Solids 443 136 (2016).  D Vouagner, C. Coussa, C. Martinet, H. Hugueney, B. Champagnon, V. Califano, and V. Sigaev, J. Non. Cryst. Solids 353 1910 (2007).  V. N. Sigaev, S. V. Lotarev, E. V. Orlova, N. V. Golubev, V. V. Koltashev, V. G. Plotnichenko, and G. A. Komandin, Gla s. Ceram. 67 105 (2010).  R. Shaltaf, H. K. Juwhari, B. Hamad, J. Khalifeh, G. M. Rignanese, and X. Gonze, J. Appl. Phys. 115 74103 (2014).  W. Kohn and L. J. Sham, Phys. Rev. 137 (1965).  H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13 5188 (1976).  S. Okamoto, A. J. Millis, and N. A. Spaldin, Phys. Rev. Lett. 97 56802 (2006).  P. B. Ghate, Phys. Rev. 133 A1167 (1964).  Y. J. Hu, D. D. Wen, Y. Q. Jiang, Y. H. Deng, and P. Peng, Trans. Nonferrous Met. Soc. China (English Ed. 25 533 (2015).  R. E. Ryltsev, B. A. Klumov, N. M. Chtchelkatchev, and K. Y. Shunyaev, J. Chem. Phys. 145 (2016).  A. Foroughi, H. Ashuri, R. Tavakoli, D. Sopu, M. Stoica, and J. Eckert, (2017).  G. Duan, M. L. Lind, M. D. Demetriou, W. L. Johnson, W. A. Goddard, T. agin, and K. Samwer, Appl. Phys. Lett. 89 23 (2006).  T. L. Cheung and C. H. Shek, J. Alloys Compd. 434 435 71 (2007).  G. Kumar, T. Ohkubo, T. Mukai, and K. Hono, Scr. Mat er. 57 173 (2007).  G. Duan, K. De Blauwe, M. L. Lind, J. P. Schramm, and W. L. Johnson, 58 159 (2008).
172  N. Barekar, P. Gargarella, K. Song, S. Pauly, U. Khn, and J. Eckert, J. Mater. Res. 26 1702 (2011).  T. Fujita, P. F. Guan, H. W. Sh eng, A. Inoue, T. Sakurai, and M. W. Chen, Phys. Rev. B Condens. Matter Mater. Phys. 81 1 (2010).  A. R. Yavari, Nature 439 405 (2006).  H. W. Sheng, W. K. Luo, F. M. Alamgir, J. M. Bai, and E. Ma, Nature 439 419 (2006).  A. L. Greer, M ater. Today 12 14 (2009).  F. Li, X. J. Liu, and Z. P. Lu, 85 147 (2014).  C. Tang and C. H. Wong, Intermetallics 58 50 (2015).  J. C. Ye, J. Lu, C. T. Liu, Q. Wang, and Y. Yang, Nat. Mater. 9 619 (2010).  Y. Shi and M. L. Falk, Phy s. Rev. Lett. 95 1 (2005).  A. J. Cao, Y. Q. Cheng, and E. Ma, Acta Mater. 57 5146 (2009).  Y. Cheng and E. Ma, Phys. Rev. B 80 064104:1 (2009).  Y. Shi and M. L. Falk, Scr. Mater. 54 381 (2006).
173 BIOGRAPHICAL SKETCH Brian Demaske was born and raised in the suburbs of Cleveland, Ohio. He moved to Florida with his family in the s ummer of 2005. He attended the University of Florida (UF) after graduating high school, but decided to transfer to the University of South Flor ida (USF) in the s p ring of 2008. In 2010, he gra duated from USF with a B.S. in p hysics. He sta rted his graduate education in p hysics the following year. There he was introduced to the field of computational m aterials s cience. His main projects include d mol ecular dynamics simulations of laser matter interactions and shock compression of crystalline metals. After graduating from USF with his M.S. in p hysics in August 2012, he applied and was accepted into the doctoral program in M aterials Science and Engineer ing at UF under the direction of Prof. Simon Phillpot co leader of the Florida Laboratory for Advanced Materials Engineering Simulation. During his graduate studies, Brian performed research in several different areas ranging from defect and diffusion beh avior in carbides to high strain rate mechanical behavior of amorphous metals to ferroelectri city in ceramics. In December 2017, he graduated from UF with a Ph.D. in m aterials s cience and engineering