USFDC Home  USF Electronic Theses and Dissertations   RSS 
Material Information
Subjects
Notes
Record Information

Full Text 
PAGE 1 Density Functional Theory Studi es of Energetic Materials by Michael W. Conroy A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy Department of Physics College of Arts and Sciences University of South Florida Major Professor: Ivan I. Oleynik, Ph.D. Venkat R. Bhethanabotla, Ph.D. Dale E. Johnson, Ph.D. Lilia M. Woods, Ph.D. Date of Approval: September 17, 2009 Keywords: equation of state, first principl es, shear stress, explosives, compression Copyright 2009, Michael W. Conroy PAGE 2 For Kathy, Larry, Corey, and Colleen PAGE 3 ACKNOWLEDGMENTS I would like to thank my research advisor, Dr. Ivan I. Oleynik, for his guidance and support during my graduate studies. He has provided me with several excellent learning opportunities and enc ouraged me towards achievements well beyond my initial aspirations as a graduate student. I thank my graduate committee, Dr. Venka t Bhethanabotla, Dr. Dale Johnson, and Dr. Lilia Woods, and the collaborators for the research, Dr. Carter T. White and Dr. Sergey V. Zybin, for their help and guidan ce. Special thanks to Mikalai Budzevich, Aaron Landerville, and Dr. You Lin for their contributions to the research and their assistance with the dissertation. I am grateful for the support from the Fred L. and Helen M. Tharp Fellowships awarded through the Department of Physics at the USF. I also thank the American Society of Engineering Education and the Na val Research Laboratory for the research experience through the NREIP. The research projects were funded by the Office of Naval Research, Naval Research Laboratory, the Army Research O ffice MultiUniversity Research Initiative on Insensitive Munitions, and Teragrid computational facilities. PAGE 4 i TABLE OF CONTENTS LIST OF TABLES ............................................................................................................. ivLIST OF FIGURES .............................................................................................................vABSTRACT ..................................................................................................................... ixCHAPTER 1 INTRODUCTION ........................................................................................11.1 Motivation and Research Objective ............................................................... 11.2Energetic Materials ........................................................................................ 51.2.1PETN............................................................................................... 61.2.2HMX ............................................................................................... 71.2.3RDX ................................................................................................ 81.2.4Nitromethane................................................................................... 91.2.5NEST1 ......................................................................................... 101.3Sensitivity to Shock ..................................................................................... 111.4Theory in EM Research ............................................................................... 121.5Outline.......................................................................................................... 13CHAPTER 2 FIRSTPRINCIPLES CALCULATIONS ..................................................142.1DensityFunctional Theory (DFT) ............................................................... 15 PAGE 5 ii 2.1.1Changing the Approach to the ManyBody Problem ................... 162.1.2Functionals of Exchange and Correlation ..................................... 182.2Methods for Practical Calculations .............................................................. 192.2.1SelfConsistent Solution of KS Equations .................................... 202.2.2Wavefunctions of Independent Electrons ..................................... 212.2.3Calculation of Physical Quantities ................................................ 24CHAPTER 3 VAN DER WAALS INTERACTIONS AND DFT ...................................263.1Dispersion Interactions ................................................................................ 273.2Empirical van der Waals Correction to DFT Calculations .......................... 28CHAPTER 4 HYDROSTATICCOMPRESSION SIMULATIONS OF ENERGETIC MATERIALS .........................................................................................354.1Previous Work ............................................................................................. 364.1.1Experimental ................................................................................. 364.1.2Theoretical .................................................................................... 384.2Computational Details ................................................................................. 404.3Equilibrium Properties ................................................................................. 424.4HydrostaticCompression Results ................................................................ 434.4.1Isothermal Equation of State ......................................................... 444.4.2Lattice Changes Under Pressure ................................................... 494.5Discussion .................................................................................................... 57 PAGE 6 iii CHAPTER 5 UNIAXIALCOMPRESSION SIMU LATIONS OF ENERGETIC MATERIALS ................................................................................................................615.1Computational Details ................................................................................. 655.2Stressstrain relationships ............................................................................ 665.2.1Principal stresses ........................................................................... 665.2.2Shear stresses ................................................................................ 815.3Band Gaps .................................................................................................... 995.4Discussion .................................................................................................. 106CHAPTER 6 CONCLUSION ........................................................................................108REFERENCES ................................................................................................................112ABOUT THE AUTHOR ....................................................................................... End Page PAGE 7 iv LIST OF TABLES Table I: Lattice Constants of PETNI, HMX, RDX, Nitromethane, and NEST1 .........................................................................................................11Table II: Homoatomic parameters of the empirical vdW correction. ......................... 34Table III: Comparison of Experimental and Calculated Equilibrium UnitCell Volumes ....................................................................................................... 43Table IV: Bulk moduli of PETNI, HMX, RDX, Nitromethane, and NEST1. .... 56 PAGE 8 v LIST OF FIGURES Figure 1: Unit cell of PETNI. .............................................................................................6Figure 2: Unit cell of HMX. .............................................................................................7Figure 3: Unit cell of RDX. ..............................................................................................8Figure 4: Unit cell of solid nitromethane. ............................................................................9Figure 5: Unit cell of NEST1. ..........................................................................................10Figure 6: The empirical vd W pair potential for interacting nitrogen atoms. .....................31Figure 7: Isothermal EOS for PETNI ...............................................................................44Figure 8: Isothermal EOS for HMX. ..............................................................................45Figure 9: Isothermal EOS of RDX. ................................................................................46Figure 10: Isothermal EOS of solid nitromethane. ............................................................47Figure 11: Isothermal EOS of NEST1. .............................................................................48Figure 12: Lattice constants of P ETN under hydrostatic compression. .............................49Figure 13: Lattice parameters of HMX under hydrostatic compression. ...........................50Figure 14: Lattice parameters of RDX under hydrostatic compression. ...........................51Figure 15: Lattice constants of solid n itromethane as functions of hydrostatic pressure. ...........................................................................................................52Figure 16: Calculated lattice constants a and c of NEST1 as a function of pressure. ......53Figure 17: Calculated lattice constant b of NEST1 as a function of pressure. .................54 PAGE 9 vi Figure 18: Principal stress as a function of volume ratio for seven lowindex uniaxial compressions ......................................................................................67Figure 19: Principal stress as a function of volume ratio. .............................................68Figure 20: Principal stress as a function of volume ratio for uniaxial compression in seven lowindex crystallogr aphic directions of PETN. ...............................69Figure 21: Principal stress (longitudinal stress) of HM X as a function of volume ratio for uniaxial compression in seven lowindex crystallographic directions. .........................................................................................................70Figure 22: Principal stress of HMX as a function of volume ratio for uniaxial compression in seven lowindex crystallographic directions. .........................71Figure 23: Principal stress of HMX as a function of volume ratio for uniaxial compression in seven lowindex crystallographic directions. .........................72Figure 24: Principal stress (longitudinal stress) of RDX as a function of volume ratio for uniaxial compression in seven lowindex crystallographic directions. .........................................................................................................73Figure 25: Principal stress of RDX as a function of volume ratio for uniaxial compression in seven lowindex crystallographic directions. .........................74Figure 26: Principal stress of RDX as a function of volume ratio for uniaxial compression in seven lowindex crystallographic directions. .........................75Figure 27: Principal stress (longitudinal stress) of nitromethane as a function of volume ratio for uniaxial compression in seven lowindex crystallographic directions. .........................................................................................................76 PAGE 10 vii Figure 28: Principal stress of nitromethane as a function of volume ratio for uniaxial compression in seven lowi ndex crystallographic directions. ...........77Figure 29: Principal stress of nitromethane as a function of volume ratio for uniaxial compression in seven lowi ndex crystallographic directions. ...........78Figure 30: Calculated maximum principal stress of NEST1 using vdWDFT for uniaxial compressions. .....................................................................................79Figure 31: Calculated maximum principal stress of NEST1 using vdWDFT for uniaxial compressions. .....................................................................................80Figure 32: Calculated maximum principal stress of NEST1 using vdWDFT for uniaxial compressions. .....................................................................................81Figure 33: Maximum shear stress for uniaxial compressions of PETN. .......................83Figure 34: Maximum shear stress for uniaxial compressions of PETN. .......................84Figure 35: Maximum shear stress for uniaxial compressions of PETN. .......................85Figure 36: Shear stress from uniaxialcompression simulations for directions of known sensitivity in PETN. .............................................................................86Figure 37: Shear stress from uniaxialcompression simulations for directions of known sensitivity in PETN. .............................................................................87Figure 38: Maximum shear stress from uniaxialcompression simulations of HMX. ...88Figure 39: Maximum shear stress from uniaxialcompression simulations of HMX. ..89Figure 40: Maximum shear stress from uniaxialcompression simulations of HMX. ..90Figure 41: Maximum shear stress from uniaxialcompression simulations of RDX. ....91Figure 42: Maximum shear stress from uniaxialcompression simulations of RDX. ...92Figure 43: Maximum shear stress from uniaxialcompression simulations of RDX. ...93 PAGE 11 viii Figure 44: Maximum shear stress from uniaxialcompression simulations of nitromethane. ...................................................................................................94Figure 45: Maximum shear stress from uniaxialcompression simulations of nitromethane. ...................................................................................................95Figure 46: Maximum shear stress from uniaxialcompression simulations of nitromethane ....................................................................................................96Figure 47: Maximum shear stress for uniaxial compressions of NEST1 using vdWDFT. ........................................................................................................97Figure 48: Maximum shear stress for uniaxial compressions of NEST1 using vdWDFT. ........................................................................................................98Figure 49: Maximum shear stress for uniaxial compressions of NEST1 using vdWDFT. ........................................................................................................99Figure 50: Band gap of PETN from uniaxialcompression simulations. .........................100Figure 51: Band gap of HMX from uniaxialcompression simulations. .........................102Figure 52: Band gap of RDX from un iaxialcompression simulations. ..........................103Figure 53: Band gap of nitromethane fr om uniaxialcompression simulations. ..............104Figure 54: Band gap of NEST1 from uniaxialcompression simulations. .....................105 PAGE 12 ix DENSITY FUNCTIONAL THEORY STUD IES OF ENERGETIC MATERIALS MICHAEL W. CONROY ABSTRACT Firstprinciples calculations employing density functional theory (DFT) were performed on the energetic materials PETN, HMX, RDX, nitromethane, and a recently discovered material, nitrate es ter 1 (NEST1). The aims of the study were to accurately predict the isothermal equation of state for each material, improve the description of these molecular crystals in DFT by introducing a correction for dispersi on interactions, and perform uniaxial compressions to investigate physical propertie s that might contribute to anisotropic sensitivity. For each system, hydrostaticcompression simulations were performed. Important properties calculated from the simulations such as the equilibrium structure, isothermal equation of state, and bulk moduli were comp ared with available experimental data to assess the agreement of the cal culation method. The largest co ntribution to the error was believed to be caused by a poor description of van der Waals (vdW) interactions within the DFT formalism. An empirical van der Waals correction to DFT was added to VASP to increase agreement with experiment. The average agre ement of the calculated unitcell volumes for six energetic crystals improved from a pproximately 9% to 2%, and the isothermal PAGE 13 x EOS showed improvement for PETN, HMX, RDX, and nitromethane. A comparison was made between DFT results with and without the vdW correction to identify possible advantages and limitations. Uniaxial compressions perpendicular to seven lowindex crystallographic planes were performed on PETN, HMX, RDX, n itromethane, and NEST1. The principal stresses, shear stresses, and band gaps for each direction were compared with available experimental information on shockinduced se nsitivity to determine possible correlations between physical properties a nd sensitivity. The results for PETN, the only system for which the anisotropic sensitivity has been thoroughly investigated by experiment, indicated a possible correlation between maximum shear stress and sensitivity. The uniaxial compressions that corresponded to th e greatest maximum shear stresses in HMX, RDX, solid nitromethane, and NEST1 were id entified and predicted as directions with possibly greater sensitivity. E xperimental data is anticipat ed for comparison with the predictions. PAGE 14 1 [Manual counter reset: Chapter 1, Section 0, Equation 0] CHAPTER 1 INTRODUCTION 1.1 Motivation and Research Objective The United States Department of Defense is actively pursuing the development of insensitive munitions to reduce the risk for loss of life and property by accidental explosions of EMs. While methods for safe handling and transportation of EMs applications are effectively utilized, the re duction of sensitivity would greatly diminish the cost of such methods. An obstacle to the de sign of insensitive materials is the fact that the microscopic events that lead to de tonation are not well understood. Also poorly understood are the mechanisms responsible for anisotropic sensitivity of EMs such as PETN, which, if clarified, mi ght aid in grasping the mechanisms of initiation. While modeling and simulation offer the ability to ca pture the miniscule time and length scales associated with initiation, modern tec hniques exhibit disadvantages ranging from prohibitive computational expense to inaccurate descriptions of atomiclevel interactions. This study aims to address three importa nt problems in EMs research: (1) firstprinciples calculations of isot hermal equations of state for important energetic materials, (2) an improvement to the description of inte rmolecular interactions in density functional theory calculations of energetic solids in the form of molecular crystals, and (3) PAGE 15 2 simulations of uniaxial compression aimed at understanding anisotropic sensitivity in singlecrystal energetic materials. The first focus of this work is to calcula te from firstprinciples the isothermal equations of state (EOSs) for important ener getic materials. The EOS of an energetic material is a very important quantity obtai ned from hydrostatic compression, allowing the calculation of the bulk modulus and its pressure derivative. From these quantities, the accessible states of a material, also termed th e Hugoniot locus, after a weak shock can be calculated [1]. Accurate and reliable firstprinciples predictions of the EOS for energetic materials are highly desired because of both the highly reactive nature of the materials and the cost associated with experiments. The second aim of the projec t is the introduction of an empirical correction to account for dispersive van der Waals (vdW) in teractions in DFT calculations for EMs. The initial calculations of this study were performed with uncorrect ed DFT, and the unitcell volumes were calculated with undesirably large error as compared with experiment. The source of the error is beli eved to be the poor (or lackin g) description of dispersion interactions by exchangecor relation functionals commonly used in DFT calculations, which is a wellknown problem in condensedmatter theory. To correct for the missing intermolecular interactio ns, a dispersion correction is used in this study that is based on the work of Neumann and Perrin [2], which showed a significant improvement in the prediction of equilibrium unitcell volume s for molecular crystals. The method was chosen to allow for quick, efficient calculations without the sacrifice of accuracy. In this study, a DFT computer code with and without the empirical correcti on (referred to below as pure DFT and vdWDFT, respectively) is used to calculate th e equilibrium unitcell PAGE 16 3 structure and to simulate the hydrostatic compression of PETNI, HMX, RDX, solid nitromethane, and NEST1. A comparison of the results is made to evaluate the improvement in the predictive capabilities of DFT with th e empirical vdW correction. The final task is to use density functiona l theory to study the anisotropic response of these materials to compression. The uniax ial compression of PETN is an important task for the study because the anisotropic sensit ivity of this material has been studied in several experiments. By performing uniaxia l compressions in the known sensitive and insensitive directions of P ETN, it might be possible to gain an understanding of the physical properties that contri bute to the anisotropy. Similar simulations and analysis are applied to HMX, RDX, nitromethane, and NE ST1, important EMs that have not been investigated for anisotropic se nsitivity to shock, to predic t the anisotropic response to compression. The anisotropic c onstitutive relations hips between stress and strain for these materials calculated from the simulations can also be useful for simulations of EM grains at greater length scales. Further, DFT is a suitable tool for this study because it has shown success in the description of EMs at high pressure at a r easonable computational expense. An important result of the simulations in this work is the prediction of the EOS for NEST1, a recently discovered EM for which hydrostaticcompression experiments have yet to be performed. Further, the EO Ss for several important EMs have been calculated in exceptional agreement with experi mental results, demonstrating the utility of DFT calculations in the pred ictive characterization of EMs. PAGE 17 4 Another significant outcome of this wo rk is that the empirical dispersion correction reduced the average error in the prediction of equilibrium unitcell volumes from 9% by pure DFT to about 2%. In addi tion, vdWDFT provides, on average, better agreement with experiment than pure DFT in the prediction of the isothermal equation of state, as well as the lattice constants as functions of pressure. These findings support the use of vdWDFT as an effective tool for use in EMs research. An additional major finding of this work is the possible existence of a correlation between shear stress and sensitivity in PETN. From the uniaxialcompression simulations, the sensitive directions, <110> and <001>, showed gr eater shear stressvalues than the insensitive directions, < 100> and <101>, at high compression. Further, the insensitive directions exhibited nonm onotonic dependence of shear stress on strain. For the other EMs of the study, compression direc tions that exhibited behavior similar to the sensitive and insensitive directions of P ETN were identified. If the sensitivity of these materials is linked with the shear stresses in a similar manner as obs erved in PETN, this work could serve as a prediction for direc tions of greater (and/or) lesser sensitivity. In the following sections, the energetic ma terials used in this study are introduced along with their properties that are relevant for this study. Th e utility of theoretical tools for the study of EMs is also briefly noted. PAGE 18 5 1.2 Energetic Materials One of the key properties of materials employed in chemical explosives is the ability to detonate. Detonation is a combusti on process that travels in the form of a shockwave at a supersonic speed through an explosive material. The chemical reactions behind the selfsustaining detona tion front release gaseous pr oducts at high pressure and temperature, typically resu lting in destructive effect s on the surroundings. Hence, energetic materials capable of detonation ar e utilized for the purposes of commercial demolition as well as warfare. Energetic materials (EMs) store large am ounts of chemical energy that can be released. The EMs used in commercial and military explosives applications may be classified as primary and s econdary explosives. Primary explosives detonate when subject to relatively weak stimuli in the form of heat or shock. A secondary explosive is typically more powerful than a primary explos ive, but will require much stronger stimuli for initiation of detonation. In a typical appl ication, a primary explos ive will be used to initiate detonation in a more powerful secondary explosive. Secondary explosives are typically form ed by organic molecules. Under ambient conditions, these materials condense to form mo lecular solids or liquids. Of particular interest for this work are the molecular cr ystals formed by the secondary explosives PETN, HMX, RDX, nitromethane, and NEST1. PAGE 19 6 1.2.1 PETN Pentaerythritol tetranitrate (PETN) is an explosive ni tric ester with chemical formula C5H8N4O12 and is known to form a molecular crystal with three polymorphs. The polymorph that is stable at ambient conditions is PETNI, and a phase transition [3] to the second polymorph PETNII occurs at a temperat ure [4] of 403 K. The discovery of a third polymorph PETNIII was recently made at high pressure [5]. PETNI, the stable polymorph at ambient conditions, forms with space group P 21c and the lattice constants of the tetragonal unit cell [6] are provided in Table I. There are 2 molecules per unit cell, shown in Figure 1, with th e total number of atoms per unit cell equal to 58. Figure 1: Unit cell of PETNI. PAGE 20 7 1.2.2 HMX Cyclotetramethylene tetranitramine (HMX) has the chemical formula C4H8N8O8, and forms molecular solids of at least four polymorphs; , , , and HMX. The polymorph is used in explosives applications because it is both stable at ambient conditions and the least sensitive to imp act [7]. The spacegroup symmetry of the polymorph is P21/c which forms a monoclinic unit cell. The lattice constants[8] of HMX are listed in Table I. The unit cell is displayed in Figure 2, and it contains two molecules or, equivalently, 56 atoms. Figure 2: Unit cell of HMX. PAGE 21 8 1.2.3 RDX Cyclotrimethylene trinitramine (R DX) is a molecule specified by C3H6N6O6 that forms solids of three polymorphs; , , and RDX. The polymorph is stable at ambient conditions, but a phase tr ansition occurs near 4 GPa to RDX [914]. The polymorph of RDX has the Pbca spacegroup symmetry with an orthorhombic unit cell [15]. The lattice c onstants[15] of RDX are shown in Table I. The unit cell of RDX, exhibited in Figure 3, is the largest of th e study, holding 168 atoms within a total of 8 molecules. Figure 3: Unit cell of RDX. PAGE 22 9 1.2.4 Nitromethane Nitromethane molecules are relatively simp le for energetic materials, having the formula CH3NO2. Unlike the other materials of this study, nitromethane is a liquid at ambient conditions, but becomes a solid at lowe r temperature [16] and/or higher pressure [17]. Several studies [1719] indicate that solid nitromet hane does not undergo a phase transition to another polymorph for pressure s below approximately 27 GPa, but this has been debated in recent work [20, 21]. Nitromethane forms a solid with P212121 symmetry at low temperature and ambient pressure. The lattice constants of solid nitromethane are displayed in Table I. The nitromethane un it cell is shown in Figure 4, and it is the smallest of the study, holding only 28 atoms (4 molecules). Figure 4: Unit cell of solid nitromethane. PAGE 23 10 1.2.5 NEST1 In the work [22] reveal ing its discovery, Chavez et al. call the new explosive nitrate ester 1. We have shor tened the name of the material, calling it NEST1 to be concise. The NEST1 molecule has a somewhat similar shape to that of PETN, and its chemical formula is C6H8N6O16. The crystal [22] has a P21/n symmetry with a monoclinic unit cell. The lattice constants [22] of NEST1 are provided in Table I. The unit cell, illustra ted in Figure 5, is relatively large, containing 144 atoms (4 molecules). Figure 5: Unit cell of NEST1. PAGE 24 11 1.3 Sensitivity to Shock One of the common methods to detonate an explosive is through a mechanical shock wave, a very fast compression at high pressure. Energetic materials (EMs) that detonate at relatively low shock pressures are said to be more sensitive to shock, whereas EMs which require much stronger shocks to cause detonation are in sensitive to shock. Owing to the fact that explosives can be initia ted by different stimuli, such as heat, shock, friction, and electric spark, one must specify th e stimulus when discussing sensitivity. For shock compression, experiments have shown that shock sensitivity in an EM can have directional dependence [23, 24]. The energetic material PETN exhibits strongly anisotropic sensitivity to shock. Initiation occurs for compression in one crys tallographic direction at approximately 4 Table I: Lattice Constants of PETNI, HMX, RDX, Nitromethane, and NEST1. The experimental references for the struct ural data are noted in the text above. System a () b () c ( ) () () () V (3) PETNI 9.38 9.38 6.70 90.0 90.0 90.0 589.50 HMX 6.54 11.05 8.7 90.0 124.3 90.0 519.39 RDX 13.182 11.574 10.709 90.0 90.0 90.0 1633.86 Nitromethane 5.183 6.235 8.518 90.0 90.0 90.0 275.31 NEST1 8.122 23.056 8.507 90.0 113.95 90.0 1456.0 PAGE 25 12 GPa [24], while shocks in excess of 30 GPa were needed to detonate the sample in another direction [25]. Meanwhile, other dire ctions exhibit sensitivity between these extremes [24]. An understanding of the diffe rence in the physical properties of PETN when compressed in different directions might contribute to an understanding of anisotropic sensitivity in EMs. On a broade r scale, knowledge of th e contributing factors to anisotropic sensitivity might help the en ergetic materials research community gain insight into the elusive microscopi c events that precede detonation. 1.4 Theory in EM Research The use of theoretical techniques allows for the study of atomicscale events that cannot be resolved in experimental measur ements. Further, the freedom of simulation design allows researchers to manipulate syst ems in ways that are difficult, if not impossible, to accomplish via experiment. Clear advantages of theory also include the reduced expense of research equipment and the use of virtual explos ives without risk of accidental explosion. As discussed above, the response of EMs to mechanical compression is an active area of EMs research. A quantummechanical th eoretical tool, density functional theory (DFT), has been successful in the descri ption of EMs at extreme conditions with exceptional accuracy, and the sp eed of modern computational technology has allowed this success to be achieved at a reasonable computational expense. On the other hand, the description of explosive molecular crystals by DFT at ambient conditions does not show PAGE 26 13 the same quality in agreement with experime ntal data. It is comm only believed that the deficiency of DFT for these materials is the poor description of weak intermolecular interactions. An ambition in modern DFT resear ch is to improve the predictive ability for systems with weak interactions, where th e successful accomplishment of the goal will have a great impact on EMs research, along w ith many other research fields involving soft matter. 1.5 Outline In Chapter 2, an overview is provided of practical densityfunctionaltheory calculations, covering the salien t features of the methods employed in this work. In Chapter 3, a description is given of the me thod used to empirically account for van der Waals interactions in DFT calculations. Chapter 4 provides an account of hydrostaticcompression simulations, as well as a compar ison of DFT calculations with and without the empirical vdW correction for each EM of the study. Chapter 5 discusses the uniaxialcompression simulations and the calculated an isotropic properties that might correlate with sensitivity to shock. In Chapter 6, the conclusions of the research are summarized. PAGE 27 14 [Manual chapter break] CHAPTER 2 FIRSTPRINCIPLES CALCULATIONS Firstprinciples calculations approach the study of systems at the atomic level from the fundamental equations of quantum me chanics rather than fitting parameters to empirically describe the system. Hence, firs tprinciples calculations are also called ab initio methods, meaning literally from the beginn ing. The most widely used theoretical tool implemented in firstprinciples calcula tions for the ground state of condensedphase systems is densityfunctional th eory (DFT). Approaches based on DFT have been used in several areas of res earch, including solids tate and biological physics, chemistry, and materials science. DFT is an attractive tool because it enables the treatment of the manybody electronic problem with feasible co st and reasonable acc uracy despite the approximations made in the handling of el ectronelectron interactions. Further, the advance of computational technology has al lowed firstprinciples DFT calculations to accurately describe systems composed of hundreds of atoms. The method was such a scientific breakthrough that the 1998 Nobel Prize in Chemistry was awarded to Walter Kohn, who provided the recipe for pr actical calculations using DFT. In this study, firstprinciples DFT cal culations were perf ormed to provide a quantummechanical description of energetic materials (EMs) at extreme conditions. This chapter is aimed to provide the reader with a concise picture of th e ideas involved in the theoretical framework of DFT and the practical methods for firstprin ciples calculations PAGE 28 15 with a specific emphasis on the methods em ployed in the project. Readers with an interest in a thorough discu ssion of the methods described below are directed to references with more information on as pects of the firstprinciples methodology. 2.1 DensityFunctional Theory (DFT) The study of a system at the atomic level begins with the consideration of a system of interacting electrons and nucl ei. Contributions to the energy of this fundamental system include the kinetic energy of the electrons and nuclei in addition to the potential en ergy describing the nucleinuc lei, electronnuclei, and electronelectron interactions ( , and respectively). The Schrodinger equation for this system is (2.1) where the manybody wavefunction is an eigenstate of the Hamiltonian operator that yields the energy of the system as an eigenvalue. Following the BornOppenheimer approximation [26], one neglects the kinetic energy of the relatively massive nuclei. Consequently, the electrons are in the instantaneous ground state corresponding to the geomet rical arrangement of the surrounding nuclei, and the electrons interact with the static external potential of the nuclear configuration. With this consideration, the electronnuclei in teraction can be e xpressed as [27] PAGE 29 16 (2.2) In addition, the Coulomb in teraction between nuclei can in many cases be treated classically, and corresponding energy term, referred to as below, is calculated with relative ease. The electronic terms require a QM treatment and an approach to the solution is to find the manybody wavefunction for the system of interacting electr ons in an external potential. A monstrous difficulty with this problem is that the wavefunction depends on the position of each electron, where the num ber of coordinates in the wavefunction depends exponentially on the number of electrons in the system. This approach is simply not feasible for systems composed of hundreds of atoms. 2.1.1 Changing the Approach to the ManyBody Problem Densityfunctional theory avoids the troubles of the manybody wavefunction by instead approaching the problem with the electron density of the system, a function of only a single pos ition. The elimination of the dependence on the coordinates of each electron and the potential for linear de pendence of calculation cost on system size for DFT provides a tremendous advantage over the manybody wavefunction approach [28]. Hohenberg and Kohn found that the gr oundstate energy of a system of interacting electrons in an external potentia l is a unique functional of the electron density PAGE 30 17 [29]. A functional maps a function to a number, and in th is case the positiondependent function representing the electron density is mapped to the energy via an energy functional (2.3) One can find the exact groundstate energy and the corresponding electron density by minimizing the energy functional with respect to the electron density (taking care to preserve the tota l number of electrons in the minimization) [28] (2.4) The missing but vital information is the mathematical form of the functional to find the exact groundstate energy. Combining aspects of the work of Hohenberg and Kohn with the work of Kohn and Sham, the functional can be expressed as (2.5) The first term represents the kinetic energy of independent el ectrons (with a factor of 2 for spin degeneracy), the second term is the classical Coulomb intera ction of the electron density (also termed the Hartree energy), and the third term in the electronnuclei interaction discussed above. The fourth term is the exchangec orrelation (xc) term that is defined are whatever terms are necessary to make (2.5) equal to the exact groundstate energy. Hidden in are the terms that are difficult and/or unknown. These include the effects of both exchange and correlation, which contribute to both the kinetic and PAGE 31 18 potential energies of the electrons [28]. Becau se the mathematical form of the exact xc functional is not known, approximations for this functional are necessary to make use of the method. Fortunately, these approximations are still able to predict groundstate materials properties with reasonable accuracy. An important step to lead to practical calculations using this theoretical framework was developed by Kohn and Sham [30] In their work, the interacting system is replaced by a system of noninteracting elec trons, as observed in (2.5), with the same density as the real system. The manybody prob lem is therefore reduced to singleparticle wavefunctions, where accounts for the missing manybody effects. 2.1.2 Functionals of Exchange and Correlation The initial xc functional used in DFT calculations was the local density approximation (LDA), where the exchange and correlation of a system is approximated by that of a homogenous electron gas. For an electron density without rapid variation, the xc functional of the LDA can be expressed as an integral [31] (2.6) where represents the energy per electron due to exchange and correlation. Note that the xc energy density, is a sum of the exchange and correlation energies. The LDA was later improved to account for spin in the local spindensity approximation (LSDA), making the xc energy density a f unction of the density for each spin, PAGE 32 19 Going a step beyond using the local approxi mation, corrections to the LSDA functional for both exchange and correlation have been developed to include the magnitude of the gradient of the density With these corrections, the xc energy density takes the form Several gradientcorrected f unctionals exist, each showing strengths for certain types of systems (s ee Refs. [27, 31] and references therein) The explicit mathematical expressions used for exch ange and correlation, as found in Ref. [27] (Chapter 8 and Appendix B), are quite complex. For the purposes of this study, it suffices to mention that the expr essions for exchange and correlation depend on the local density, where the addition of th e density gradient dependence introduces a semilocal approximation. Meanwhile, the weak intermolecular interact ions prevalent in molecular crystals are insufficiently descri bed by these functionals because of their nonlocal character. The method introduced to remedy the DFT calculations in this study is discussed in Chapter 3. 2.2 Methods for Practical Calculations To make use of DFT, approximations and other tricks must be introduced to make the solution process feasib le. Several of these important considerations for practical DFT calculations are briefly discussed belo w. For a more thorough discussion of the methods used in the solution of the KS equations, see Ref [27]. PAGE 33 20 2.2.1 SelfConsistent Solution of KS Equations The independentelectron equations of Kohn and Sham must be solved numerically by finding the electr on density and effective potentia l that are consistent. In the solution process [27], an initial guess is made for the spindependent density of electrons, and From this input density, an e ffective potential is calculated from the electron density [27], (2.7) and used in the KohnSham (KS) equation [27], (2.8) to solve for independentelectron states. Note that the dependence of the potential terms on spin have been denoted with the superscript From the independentelectron wavefunctions, the output electron density may be calculated as (2.9) where denotes the occupation of the correspond ing state [27]. At this step in the solution process, the guessed input density ma y be compared with the output density to check for selfconsistency. In the event that input and output dens ities are not equal within a desired tolerance, a strategic guess (see Ref. [27], Chapter 9 for details) is made for the subsequent input density. When the input density su fficiently matches the PAGE 34 21 output density, selfconsistency is reached. Fr om the selfconsisten t solutions, the output quantities of interest such as energy, forces stress, and eigenval ues are calculated. Reaching selfconsistency does not truly occu r in practical calcul ations because it is both expensive and unnecessary to conti nue the solution process until the input and output densities are exactly the same. For practical purpos es, the densities must match within a tolerance value, usually specified by a user within DFT computer codes. This tolerance is typically in the form of an energy value, wher ein the solution is acceptable when the difference in energy between the i nput and output densitie s is less than the specified tolerance. 2.2.2 Wavefunctions of Independent Electrons DFT requires the considerat ion of the electron density rather than the manybody wavefunction, but wavefunctions are not disregarded comple tely. The calculation of independentelectron wavefuncti ons is clearly needed for th e solution of the KohnSham (KS) equations. Different approaches exist for the solution of the KS equations, and a manner of classifying the methods is through th e type of basis used for expansion of the independentelectron wavefuncti ons. The expansion of a wavefunction in a basis can be performed in a number of ways, but a convenient choice for solids is to have an orthonormal basis set that obeys periodic bounda ry conditions. Plane waves are a widely used basis set for DFT calculations of solids, and their use provides several advantages. Some of the benefits include completeness, the same treatment of all space, the PAGE 35 22 independence of the basis on atomic positions, simplicity in mathematical form and in the calculation of derivatives in re ciprocal space, and the ease by which convergence can be achieved [32]. The Fourier expansion of such a wavef unction yields the desired planewave expansion [27], (2.10) Exploiting the orthonormality of the basis se t yields a Schrdinger equation for each eigenstate [27], (2.11) By defining and where is a lattice vect or in reciprocal space, the Schrdinger equation can be written for a point in reciprocal space as [27] (2.12) A difficulty that must be faced is the large number of plane waves that are required to describe a rapidly varying wavefunction. Such vari ations occur in the region of the core electrons, near the nucleus of an atom. Fortunately, pseudopotential methods have been developed, and a benefit of thes e methods is that the variation of the wavefunction is smoothed out near the nuc leus. As a result, fewer plane waves are required for the calculations. The pseudopotenti als smooth the variations by forming an effective potential that incorporates the eff ects of both the core electrons and nuclei. PAGE 36 23 Because of the relatively insignificant contribut ion of the core states to chemical bonding and other properties of the solid state [33] pseudopotentials are great approximations for DFT calculations that possess the additional advantage of reducing the computational expense of calculations. A number of met hods exist to form these potentials, and a summary of these techniques may be found in Ref. [27], Chapter 11. With plane waves, the convergence of a physical quantity calculated from the DFT formalism can be achieved simply by increasing the number of plane waves in the basis until the quantity no longer varies or, in practice, vari es less than a desired error with increased cutoff. With the proper treatment to produce a wavefunction without rapid variation, the wavefunction can effectively be expanded in a finite number of plane waves. Noting that regions of rapid varia tion of the wavefunction are characterized by high kinetic energy, a maximum kine tic energy, typically called the cutoff energy, can be set to indicate the magnitude of the maximu m wavevector used in the expansion of the wavefunction. The cutoff energy and the maximum wavevector are related by the expression [33] (2.13) which is the energy of a free electr on with the corresponding wavevector. PAGE 37 24 2.2.3 Calculation of Physical Quantities The calculation of several important physic al quantities such as the energy per unit cell requires the integrati on over a continuous variable within the Brillouin zone. Fortunately, these integrals can effectively be reduced in most cases to an average over a small number of points without much loss on accur acy. The number of points needed for accurate integration is great er for integrands with rapid variation. For insulators such as those of this study, the ba nds are filled, and integrati on can be done with a small number of special points. For metals, ma ny more points are needed to sufficiently sample the Fermi surface, where the occupation of bands changes quickly. There are several methods for the calcul ation of special points for the use of integration within the Brilloui n zone [34]. Overall, the met hods are aimed to use crystal symmetry to approximately find the point where the integral and inte grand are equal [27]. A popular choice is the MonkhorstPack method [35], which is convenient because the special points can be found for any crysta l type with a simple formula [27]: (2.14) The sum is performed over three dimensions, where are the lattice vectors in reciprocal space. These points form a uniform grid of points in reciprocal space that is translated off of the point. The forces and stress in a system are im portant physical quantities to extract from DFT calculations. The force conjugate to a ny parameter within the Hamiltonian of a PAGE 38 25 system may be calculated with the HellmanFeynman theorem [27, 36]. The force on an atom and the stress tensor from firstprinci ples can be found from derivatives of the energy with respect to atomic coordinates and strain components, respectively. The force on atom can be found by the derivative of the en ergy with respect to the coordinates of atom [27]: (2.15) Separating the terms in the energy that requ ire a quantummechanical treatment into the Hamiltonian and the energy of the classical nucleinuclei interaction the expression for the force on an atom in the ground state can be written as [27] (2.16) This expression is a result of firstorder perturbation theor y, where the derivatives of the wavefunction with respect to th e atomic coordinate vanish in the groundstate because the wavefunction is at extrema with respect to variations in the parameters describing the system [27]. Similarly, a sy stem in equilibrium will have a stress calculated by taking the derivative of the ener gy with respect to strain per unit volume, (2.17) where is given by the scaling to The numerical differentiation of the energy to calculate the forces and stress in this study is accomplished within the rou tines of the VASP computer code [37]. PAGE 39 26 [Manual chapter break] CHAPTER 3 VAN DER WAALS INTERACTIONS AND DFT The commonly used exchangecorrelati on functionals of the local density approximation (LDA) and the generalized gr adient approximation (GGA) do not yield the correct asymptotic behavior of longrange interactions. The improper description of these interactions by DFT gives poor results fo r the structure of sparse systems, where longrange interactions are vital for structur e determination. Examples of sparse systems that are poorly described by DFT include mol ecular crystals, vdW complexes, proteins, DNA, and other biomolecules. Energetic materi als, which are usually molecular crystals at ambient conditions, fall into this category, and DFT results for these systems with LDA and GGA functionals are insufficiently accurate [38]. For example, the unitcell volumes of energetic materials are predicted with an error of roughly 10% [38]. Owing to the fact that several important properties of materials such as the detonation velocity are functions of the density (or powers of the de nsity), the error in vol ume prediction desired by experimentalists is less than 1% [38]. The cause of the error for sparse systems is commonly believed to be due to longrange interactions. Because the commonly used xc functionals depend on the local or semilocal density, only effects at this range are well described. For neutral systems, an important longrange interacti on that cannot be described by common DFT functionals is PAGE 40 27 dispersion, and the resulting intermolecular forces are typically called van der Waals forces or London forces. 3.1 Dispersion Interactions Dispersion forces are a quantummechanical effect where fluctuations in electron density separated in space interact with one another. These forces are a type of polarization force between molecules in that their origin is from the polarization of one molecule by another [31] and can be expr essed by secondorder RayleighSchrdinger perturbation theory (see Ref [31], Section A3.3). The dispersion interaction can be expresse d as a multipole expansion of the form [31]: (3.1) where the sum includes only even powers of for interacting atoms separated by the interatomic distance The values of are termed the dispersion coefficients. An important consideration is that the l eading term that desc ribes instantaneous dipoledipole interactions can account for most of the dispersion inte raction. In fact, at separations of about 10 bohr (~5.3 ), the er ror introduced by using the dipoledipole term is less than 20% [31], and the er ror decreases at larger separations. PAGE 41 28 3.2 Empirical van der Waals Correction to DFT Calculations Two common methods are employed to treat vdW forces in DFT calculations: the firstprinciples approach and the empirical ap proach. Firstprinciples approaches [3951] maintain the spirit of DFT calculations by avoiding any fitting to experimental data, and this approach shows great promise. However, firstprinciples methods for a dispersion correction to DFT ar e currently too expensive for large systems. Further, the results from firstprinciples a pproaches have not shown a sign ificant increase in accuracy that would justify the comput ational expense. On the othe r hand, the empirical approach [2] has recently been shown to significan tly improve the accur acy of calculations involving molecular crystals [2] at a much sm aller cost than firstprinciples methods. Although the empirical method deviates fr om the firstprinciples methodology and requires experimental data, it is currently able to predict unitcell volumes within an error that the firstprinciples methods are still striving to achieve. In this work, an empirical approach based on the work of Neumann and Perrin [2] is used for the simulation of energetic material s. For the molecular crystals in their study, the average error in the pr ediction of unitcell volume was reduced from 20.4% by pure DFT to about 1% with the em pirical vdW correction [2]. Th is result is a significant improvement over previous works on the us e of empirical vdW corrections for the prediction of equilibrium unitcell structures of molecular crystals. PAGE 42 29 Following the work of Neumann and Perrin [2], the salient points of the approach are summarized below. The total energy of th e unit cell is a sum of two parts, the energy calculated by DFT and the energy from the vdW correction [2] (3.2) The Vienna Ab Initio Simulation Package (VASP) [37, 52] is a DFT code used to calculate the energy and its derivatives for the DFT component. Although VASP is treated as a black box in the calculation method, the parameters for calculations must be chosen properly. The details of the VASP para meters used in this work are provided in Chapter 4. The dispersion component of the approach is formed by interaction potentials between pa irs of atoms [2]: (3.3) The first two sums include the atoms within the unit cell, and the third sum includes all realspace lattice vectors. The primed sum is an indication that the selfinteraction must be avoided, and the factor (1/2) is introduce d to avoid double counting. The pair potentials depend on the type of the interacting atoms and and the type depends on both the atomic number of the atom and the number of its covalent bonds. The interatomic distance is specified as the distance between atom translated by the lattice vector and atom Only atoms separated by 18 or less are considered in order to truncate the infinite sum, and a spline function brings the potentials to zero smoothly between 15 and 18 . PAGE 43 30 The pair potentials are formed by a product of two terms [2]: (3.4) The leading term in the multipole expans ion of the dispersion interaction is (referred to as the dispersion term below), and it is multiplied by a damping function that controls the distancedependent influence of the vdW pair potential At long distances, the dispersion term alone properly describe s the asymptotic behavior, and the damping function appropriately takes a value of one. At short dist ances, the damping function must prevent the divergence of the dispersion term. A breakthrough of the approach is the improvement of the damping function at intermediate distances. Because the GGA functionals PW91 [53, 54] and PB E [55] can provide an accu rate description of shortrange interactions consistent with ionic and covalent bonds, it is crucial for the vdW correction potential to properly turn on at longer distances where the GGA functional description begins to fail. The da mping function takes the form [2] (3.5) PAGE 44 31 where the form factor was introduced to add flexibility to the similar potential by Mooij et al. [56] used in a previous empirical vdW approach [57]. The parameter is the crossover distance, and it denotes the inte ratomic distance where the dispersion term intersects the constant value taken by the pair potential at very small interatomic Figure 6: The empirical vdW pa ir potential for interacting nitrogen atoms. The pair potential is shown in solid black. The da shed red line is the dispersion term and the dashed green line is th e constant value that the pair potential approaches at small interatomic distance. The intersec tion of the green and red dashed lines determines the crossover distance in the damping function. PAGE 45 32 distances. In Figure 6, the empi rical vdW pair potential for interacting nitrogen atoms is shown to illustrate the influe nce of the damping function.. The number of fitting parameters is reduced by simple rules to calculate heteroatomic parameters from homoatomic parameters. For the crossover distances, the values for different atoms of type A and B are found by the average of the homoatomic values according to [2] (3.6) The heteroatomic coefficients are calcul ated from homoatomic coefficients via the combination rule [2] (3.7) devised by Wu and Yang [57]. The effective number of electrons is represented by the parameter and these values are taken fr om the work of Halgren [58]. The empirical nature of the correction is in the fitting of the parameters , and to experimental data. Using the molecular coefficients calculated from the dipole oscillator strengths measured by Meath and co workers [5964], a form is made for the calculation of atomic coefficients assuming that molecular coefficients simply added [2]: (3.8) PAGE 46 33 The molecules are denoted by the Greek subscripts and and the molecular coefficients are formed by a sum over the coefficients for interacting atoms on the different molecules. Using this form and the combination rule, the homoatomic coefficients are found by mi nimizing the function [2] (3.9) The crossover distances and the form f actor are fit by a single quantity that accounts for the difference in both lattice c onstants and angles as compared to experimental structures [2]. The fitting procedure is covered in detail in the appendix of Ref. [2], but a brief expl anation is provided below. A transformation between pred icted and experimental latt ice vectors is cleverly formed by a product of a rotation matrix and a symmetric matrix with orthogonal eigenvectors. The eigenvalues of the symmetric matrix represent compressions along the corresponding orthogonal eigenvectors, and the rotation matrix doe s not alter the unitcell geometry. From the eigenvalues a quantity is defined as [2] (3.10) to represent the difference between the predicte d and experimental lattice. Essentially, the crossover distances and the form factor are determined by the minimization of for the molecular crystals of the study. The experimental data set incl udes the unitcell structures of 31 molecular crystals obtai ned at lowtemperature [2]. PAGE 47 34 The homoatomic parameters used for the model are shown in Table II [2, 58]. Note that heteroatomic parameters can be calculated from the values shown. The coefficients depend on the atom type, which includes both the chemical element and the number of its covalent bonds. Only one t ype was used for nitrogen because of limited availability of experimental data for mol ecules containing nitrogen [2]. The crossover distances and effective electr on numbers, however, depend only on the chemical element. For all interactions, the form factor =0.25 was used. Table II: Homoatomic parameters of the empirical vdW correction. All heteroatomic interactions can be calculated from these data and the combination rules shown above. Atom Number of covalent bonds Homoatomic coefficient (eV 6) Homoatomic crossover distance () Effective number of electrons, Carbon 2 17.87 3.884 2.49 Carbon 3 16.47 3.884 2.49 Carbon 4 13.09 3.884 2.49 Oxygen 1 7.645 2.837 3.15 Oxygen 2 6.904 2.837 3.15 Nitrogen (any) 11.41 3.384 2.82 Hydrogen 1 1.702 3.200 0.8 PAGE 48 35 [Manual chapter break] CHAPTER 4 HYDROSTATICCOMPRESSION SIMULATIONS OF ENERGETIC MATERIALS The application of hydrostatic pressure to a system is analogous to the compression of a static fluid, wherein shear stresses are absent and the compressive stresses, i.e. the diagonal st resstensor elements, are equal. The application of hydrostatic pressure to a material may be used to inve stigate its equation of state (EOS), a relation between thermodynamic state variables such as the pressure, volume and temperature of the material. The EOS of a material is of cr itical importance for a multitude of practical problems. For energetic materials, an exam ple application of the EOS from hydrostaticcompression data it allows one to calculate th e bulk modulus and its pressure derivative, which can be used to calculate the linear relationship between shock and particle velocities for weak shocks in a material [1 ]. For this study, the temperature will be restricted to 0 K because DFT is a groundst ate theory. Hence, this work involves the simulation of hydrostatic compression on PETN HMX, RDX, nitromethane, and NESTI to calculate the isothermal EOS, relating pre ssure and volume. Also investigated are the pressureinduced structural changes in the unit cell. PAGE 49 36 4.1 Previous Work The hydrostatic compressions of PETN, HMX, RDX, and solid nitromethane have been performed in both experimental and theoretical studies by other research groups. A brief review of the relevant pr evious works is provided in the following subsections. 4.1.1 Experimental Hydrostatic compression experiments on the EMs of interest have been carried out using diamond anvil cells. In the experiments, a very small sample is placed inside of a hole drilled in a metal gasket along with a pressuredistributing me dium that (ideally) creates hydrostatic conditions under pressure. The pressure is applied by diamond anvils on either side of the sample/medium mixture. Xrays can pass th rough the apparatus at desired pressures for analysis of the sample. Hydrostatic compression experiments on P ETN up to 10 GPa at room temperature were completed by Olinger et al [12]. Yoo et al. [14] also studied the hydrostatic compression of PETN, extending the pressure up to 15 GPa. Experiments of hydrostatic and nonhydr ostatic compression on HMX have been performed by Yoo and Cynn [65], and late r by Gump and Peiris [66]. Yoo and Cynn argued from their xray diffraction and Ra man spectroscopy results that HMX undergoes PAGE 50 37 a phase transition at the pre ssures of 12 and 27 GPa from HMX to and HMX, respectively [65]. A lower pre ssure regime up to 5.8 GPa was studied in the experiments of Gump and Peiris, where HMX remained in the beta phase duri ng the application of hydrostatic pressure at 30, 100, and 140 C [66] However, their samples converted to the phase upon decompression. Olinger et al. [13] and Yoo et al. [14] have reported the isothermal EOS of RDX from experiments of hydros tatic compression. Both studies observed phase transitions near 4 GPa from to RDX [13, 14]. The transition has been investigated with Raman spectroscopy experiments by Ciezak et al [9] and Dreger and Gupta [11]. The structure of the polymorph was recently identified by Davidson et al [10]. The EOS has been reported from experi ments of hydrostatic compression on solid nitromethane done by Cromer et al. [17], Yarger and Olinger [19], and Citroni et al. [18]. Cromer studied compression up to 6.0 GPa, Yarger and Olinger up to 15 GPa, and Citroni et al. up to 27.3 GPa, nearly approaching the detonation threshold pressure[1719]. Each of the studies repor t that solid nitromethane maintains its crystal symmetry throughout its respective pressure range. However, Courtecuisse et al. [20] have reported Xray diffraction results indicating solidsolid phase transitions at approximately 3, 7.5, 13.2, and 25 GPa. It has been argued [21] that these results are indi cative of changes in molecular structure rather than phase transitions. PAGE 51 38 4.1.2 Theoretical Modeling techniques have also been employed to study the hydrostatic compression of PETN, HMX, RDX, and nitromethane. DFT calculations to simulate the hydros tatic compression of PETN have been performed by Gan et al. [67] and Byrd and Rice [68]. In the work of Gan et al. the PBE xc functional was used with a Gaussian basi s set [67]. Byrd and Rice have studied the compression of PETN (as well as HMX, RDX, and other important EMs) with a planewave basis using the PW91, PBE, and LDA xc functionals at diffe rent energy cutoffs [68]. Hydrostatic compression has also been modeled at the HartreeFock level by Brand [69]. Several theoretical inves tigations of the hydrostati c compression of HMX have been reported in the literature. Using a rigidmolecule approximation and classical interatomic potentials fit to nearambient conditions, Sewell [70] has completed Monte Carlo simulations of beta HMX under hydros tatic pressures up to 7.5 GPa. Sorescu et al. [71] have performed molecularpacking and mo leculardynamics simulations on HMX as well as RDX, under hydrostatic compression w ith the constraint of a rigidmolecule approximation. Another study by Sewell [72] i nvestigated the compression of the pure phases of HMX via molecular dynamics calcula tions, predicting the isotherms and elastic moduli for each polymorph. HartreeFock cal culations of compression on HMX have been done by Zerilli and Kuklja [73]. DF T studies of hydrostatic compression of HMX have also been executed with the LDA xc functional up to 40 GPa [74]. PAGE 52 39 As mentioned above, Byrd and Rice [68] have performed hydrostatic compression simulations on HMX and RDX. While theoreti cal investigations of HMX are abundant, the size of the RDX unit cell is most likely a source for the low number of calculations involving hydrostatic compression for this EM. The computational expense for firstprinciples calculations of system s as large as RDX is prohibitive. In contrast to RDX, solid nitrometha ne has a very small unit cell, and is commonly regarded as a prototype EM for simulation. The hydrostatic compression of nitromethane was studied via HartreeFock ca lculations by Zerilli [75] and with DFT by Liu et al. [76] up to 20 GPa. Reed et al [77] have performed ex tensive DFT calculations on solid nitromethane, which include hydrosta tic compression up to very high pressures. The uniqueness of the calculations of th e present study for EM research is the implementation of the vdW correction into DF T codes to study these important systems both at equilibrium and under hydrostatic comp ression. The lacking pr oper description of dispersion forces in DFT has been noted fo r EMs and other molecular crystals at low pressure, and the empirical vdW correction has been shown to remedy the predictive ability of DFT for systems such as nitrometha ne [2]. However, the agreement of results with the vdW correction has not been studied for other EMs, and has not been studied for any molecular crystals under hydrostatic pr essure. Through the comparison of the vdWcorrected results with both pure DFT and experiment, the ability of the empirical correction to properly describe dispersion effects with pre ssure can be investigated. Further, the hydrostaticcompression calcula tions serve as a method to assess the predictive ability of our DFT approach because experiments have been performed with data available for comparison. However, it is not possible to co mpare the physical PAGE 53 40 quantities obtained from the uniaxialcomp ression simulations described in Chapter 5 (page 61) with experiment, simply because they either have not been performed or experimental techniques were not used (or in ex istence) to probe the quantities of interest. Some of the results below have appeared in previous publications [7880] and are shown below with permission where appropriate. 4.2 Computational Details The firstprinciples DFT [29, 30] calcu lations to model hydrostatic compression were performed for the EMs using the Vienna AbInitio Simulation Package (VASP) [37, 52]. Tests were completed to determine the pa rameters of the calculations that would provide sufficient accuracy while minimizing computational expense. Parameters that significantly contribute to the accuracy of th e results include the exchangecorrelation functional, pseudopotential, kpo int set, and energy cutoff. The combination of exchangecorrelati on functional and pseudopotential was chosen to provide the best agreement with e xperimental structure. The functionals tested were the PW91 [53, 54] and PBE [55, 81] functi onals. Both functionals were tested with the PAW potentials [82, 83] as implemented in VASP, and the PW91 functional was also tested with ultrasoft pseudopotentials (U SPs) [84, 85]. In the initial tests, each combination of functional and potential was us ed in the relaxation of the experimental structure of PETN. The parameters of each relaxation included a reasonably high energy cutoff of 1,000 eV and MonkhorstPack (MP) grid of 2x2x2 (spacing: 0.07 1). For PAGE 54 41 these calculations and all others in the study, an ener gy tolerance of 106 eV was used as the criteria for selfconsistent solution of the KohnSham equations. The PBE functional with the PAW potential showed the best agreement with experiment, and this combination was chosen for the DFT component of all calculations. Tests were subsequently performed to de termine the kpoint sampling density for the calculations. The kpoint sampling was tested first by setting the energy cutoff to the high value of 1,000 eV and executing totalenergy calculations on the experimental structure, both at 50% volume comp ression and without compression. The halfcompressed structure was constructed by scaling the latti ce constants by 0.51/3 and keeping the atoms in fractiona l coordinates. The rationale for the tests at the two compressions is that the Brillouin zone e xpands when its realspace counterpart is compressed, and the kpoint sampling density for a fixed number of kpoints will decrease with compression. The tests were ai med to be sure that the accuracy of the calculations was sufficient for all compre ssion values of the study. For PETN, it was found that an average kpoint sampling density of 0.08 A1 at half compression yielded convergence by less than 1 meV in energy per atom and 0.05 GPa in pressure. The other EMs of the study yielded similar convergen ce results, and the average sampling density of 0.08 1 at halfcompression was used to dete rmine the MP grid for the hydrostaticcompression simulations of each EM in the study. Kineticenergy cutoff tests were performed with the parameters specified above to find the minimum planewave basis set that provided sufficient accuracy. Total energy calculations of the experimental structures at ambient conditions were performed at cutoff values of 400, 500, 700, and 1,000 eV. The 700 eV calculation for PETN gave PAGE 55 42 convergence better than 0.4 GPa in pressure and 1.5 meV per atom in energy. Similar convergence values were obtained for the othe r EMs of the study at an energy cutoff of 700 eV, which was the cutoff chosen for the hydrostaticcompression simulations of the work. The chosen cutoff is 1.75 times larger than the default cutoff specified by the PAW potentials in VASP, which is us ually necessary for molecular crystals to obtain sufficient convergence. The experimental structures of each EM were relaxed via the conjugate gradient algorithm as implemented in VASP to obtain the theoretical equilibrium structure. In the energy minimization, the unit cell volume, shape, and symmetry as well as the atomic coordinates were allowed to relax until th e maximum force on any atom was less than 0.03 eV/. It is expected that the error from the choice in fo rce tolerance is much smaller than the error from the chosen energy cutoff. The hydrostaticcompression simulations we re completed by scaling the volume of the unit cell in increments of 2% from the calculated equilibrium volume V0 to V/V0 = 0.60. At each step, the conditions of hydrostatic pressure we re simulated by holding the unitcell volume constant while relaxing the at omic coordinates and unitcell shape. In the relaxations, the symmetry of the unit cell was not constrained. 4.3 Equilibrium Properties The calculated equilibrium structures fo r PETN, HMX, RDX, solid nitromethane, and NEST1 using both pure DFT and the em pirical vdW correction are displayed in PAGE 56 43 Table III along with the percent error as comp ared to the experimental volumes. The unitcell volumes predicted by pure DFT (GGAPBE) are approximately 9% larger on average than experiment. As described in Chap ter 3, it is commonly believed that a poor description of vdW forces within GGA xc func tionals is the main s ource of the error in the prediction of equilibrium structures for molecular crystals, where intermolecular forces are significant. With the vdW correction, th e calculated unitcell volumes are less than experiment by 2.3% on average. 4.4 HydrostaticCompression Results From the hydrostaticcompression simulati ons, the isothermal equation of state (EOS) and bulk moduli were calculated. The lat tice constants were al so calculated as a Table III: Comparison of Experimental a nd Calculated Equilibrium UnitCell Volumes System Experimental Pure DFT Pure DFT %error DFTvdW DFTvdW %error PETNI 589.50 (Ref. [6]) 620.11 +5.2% 567.27 3.8% HMX 519.39 (Ref. [86]) 556.07 +7.1% 500.77 3.6% RDX 1633.86 (Ref. [15]) 1775.95 +8.7% 1591.21 2.6% Nitromethane 275.31 (Ref. [16]) 304.22 +10.5% 270.00 1.9% NEST1 1456.01 (Ref. [22]) 1649.48 +13.3% 1464.18 +0.6% PAGE 57 44 function of pressure. The results from the si mulations are compared with experimental data for PETN, HMX, RDX, and nitromethane. 4.4.1 Isothermal Equation of State The isothermal EOS of PETN is show n in Figure 7 and compared with the experimental data of Olinger et al. [12]. With respect to th e experimental data, the Figure 7: Isothermal EOS for PETNI PAGE 58 45 agreement is clearly increased with the addition of the vd W correction. The EOSs from some other experiments are not provided because volume ratios (V/V0) were reported in the absence of absolute volumes and/or equilibrium volumes. The EOS for HMX is displayed in Figure 8 with the results of two hydrostatic compression experiments [65, 66]. In the lowpr essure regime, the volume of the unit cell at a given pressure is overestimated by pure DFT and underestimated by the vdW correction when compared to the experiment al data. For the study by Yoo and Cynn [65], the vdW results appear to agree well at pr essures between 525 GPa. The discrepancy Figure 8: Isothermal EOS for HMX. PAGE 59 46 between all DFT calculations a nd experiment at higher pressu res increases, most likely because of a phase transition[65] to HMX near 27 GPa. This phase transition was measured to have a 4% change in unitcell vo lume [65]. In fact, another phase transition was predicted at 12 GPa, but the resulting change in volume was negligible [65]. From the simulation data, no phase transitions were observed: the spacegroup symmetry remains the same for all structures during the compression. The EOS for RDX is provided in Figure 9, along with the data from the experiments of Olinger et al. [13]. Below the phasetransition pressure of 4 GPa, the Figure 9: Isothermal EOS of RDX. PAGE 60 47 agreement of the EOS with experiment is improved with the addition of the vdW correction. The experimental data show a reduction in unitcell volume after the transition, however the transition wa s not observed in the simulation. The EOS for solid nitromethane is displa yed in Figure 10, and experimental data from Refs. [1719] are shown for comparis on. At low pressure, the vdWcorrected DFT data appear to agree well with the ex perimental data [17, 18] of Citroni et al. and Cromer et al. but the pure DFT data show better agr eement with the data [19] of Yarger et al. At pressures above 10 GPa, however the Citroni and Yarger ex periments approach the EOS Figure 10: Isothermal EOS of solid nitromethane. PAGE 61 48 calculated with pure DFT by showing a greater pressure at a given volume. Meanwhile, the Cromer experiment shows the opposite tr end, yielding a significantly smaller pressure at a given volume. As discu ssed in Ref [19], the greater compressibility measured by Cromer et al. might be a result of the geometrical alignment of the xray beam along the axis of the diamond anvil cell. The EOS for NEST1 is shown in Figure 11. There are no known hydrostatic compression experiments for comparison. Based upon the results for the other systems, it is likely that the experimental data will reas onably agree with the results obtained from Figure 11: Isothermal EOS of NEST1. PAGE 62 49 the vdW calculations provided that NEST1 does not undergo a phase transition in the pressure interval studied. 4.4.2 Lattice Changes Under Pressure The lattice parameters of PETN as a f unction of pressure from the hydrostaticcompression simulations are compared with expe rimental data [12] in Figure 12. Below 5 Figure 12: Lattice constants of PETN under hydrostatic compression. The b lattice constant is not displayed because PETN has a tetragonal unit cell ( a = b ). PAGE 63 50 GPa, the vdWcorrected DFT data show bette r agreement with experiment than pure DFT. Above 5 GPa, the pure DFT data for the lattice constant a appears to show slightly better agreement than the vdWDFT data. On the other hand, vdWDFT clearly shows better agreement with experiment for the c lattice constant. The calculated and experimental [65, 66] lattice constants of HMX are shown in Figure 13. It is difficult from the plot to ev aluate the improvement in agreement obtained with the use of vdWDFT. The predictions for the a and c lattice constants via both methods are in similar agreement with expe riment. The vdWDFT data appear to agree Figure 13: Lattice parameters of HMX under hydrostatic compression. PAGE 64 51 slightly more with the experimental data of Yoo and Cynn, while the pure DFT look to agree more with pure DFT. A similar trend is shown for the EOS calculation in Figure 8. The calculated lattice parameters of R DX are displayed in Figure 14 along with experimental data from Ref [13]. The vdWcorrected data clearly show better agreement with experiment than pure DFT for the a and c lattice constants up to the phasetransition pressure of about 4 GPa. Meanwhile, the agreement for both calculation methods is similar for the b lattice constant up to about 4 GPa. Overall, the vdW correction yields pressuredependent lattice constants with better average agreement that pure DFT. Figure 14: Lattice parameters of RDX under hydrostatic compression. PAGE 65 52 The calculated lattice constants of solid nitromethane are displayed in Figure 15 along with experimental data [1719]. The b lattice constant calculated with the vdW correction agrees more with the experimental data than pure DFT, but the opposite is true for the c lattice constant. Meanwhile, both calculati on methods predict the lattice constant a in good agreement with experiment. Similar to the case of HMX, an evaluation of the improvement in the prediction of lattice constants of nitr omethane provided by the vdW correction is difficult to assess. The calculated lattice constants a and c of NEST1 from the hydrostaticcompression simulations are displayed in Figur e 16. Owing to the larger length of lattice Figure 15: Lattice constants of solid nitromethane as functio ns of hydrostatic pressure. PAGE 66 53 constant b it was plotted separately in Figure 17. As mentioned above, there are currently no experimental data available for comparison. Figure 16: Calculated lattice constants a and c of NEST1 as a function of pressure. PAGE 67 54 4.4.3 Bulk Moduli The bulk modulus is defined as and it is a measure of the stiffness of a material. The bulk modulus and its derivative with respect to pressure can be calculated by fitting hydrostatic compression data to equati ons of state for pressure as a function of volume. The bulk moduli for each EM of the study were calculated by fitting data to three common equations of state; the BirchMurnaghan, Murnaghan, and Vinet EOS forms. The BirchMurnaghan equation is given by Figure 17: Calculated lattice constant b of NEST1 as a function of pressure. PAGE 68 55 (4.1) the Murnaghan equation is given by (4.2) and the Vinet equation is given by (4.3) The calculated bulk moduli are shown in Ta ble IV along with reported experimental values obtained by various methods. For consis tency with reported experimental values, the hydrostaticcompression data were fit to the EOS forms up to different pressures for each system. PETN was fit up to 10 GPa, HMX up to 12 GPa, RDX up to 4 GPa, and nitromethane up to 15 GPa. The data for NEST 1 was fit up to 5 GPa. The lower pressure was chosen because the fitting forms provided very similar results in this pressure range. Also, it is not known if NEST1 experiences a phase transition to a different polymorph in the pressure range studied. By choosing a lower pressure for fitt ing, there is a better chance that the undesirable case might be avoi ded wherein the maximu m pressure chosen for fitting is above the transition pressure. PAGE 69 56 Table IV: Bulk moduli of PETNI, HMX, RDX, Nitromethane, and NEST1. The superscipts BM, M, and V indicate that th e value was obtained by fitting data to the BirchMurnaghan, Murnaghan, a nd Vinet EOS, respectively. DFT without vdW DFT with vdW Experiment System B0 B0 B0 B0 B0 B0 PETNI 9.0BM 10.4M 9.2V 8.6BM 5.9M 7.9V 14.2BM 16.0M 14.7V 9.0BM 6.4M 8.1V 12.3BM, [14] 9.9[12] 9.1[67, 87] 8.2BM, [14] 11.0[12] HMX 11.0BM 12.3M 11.2V 8.4BM 6.1M 7.9V 19.3BM 20.8M 19.5V 7.9BM 6.1M 7.6V 21.0BM, [66] 12.4BM, [14] 13.6[13] 7.45BM, [66] 10.4BM, [14] 9.3[13] RDX 10.0BM 10.4M 10.0V 7.2BM 6.2M 7.1V 15.9BM 16.4M 15.9V 7.8BM 6.6M 7.6V 13.9BM, [14], 13.0, 9.8BM 11.38[88], 11.05[88] 11.94[89], 11.67[89] 12.05[90], 11.92[90] 5.8BM, [14] 6.6[88] 11.4BM Nitromethane 6.9BM 8.0M 6.8V 7.2BM 5.1M 7.2V 11.7BM 12.9M 11.7V 7.0BM 5.3M 7.0V 8.3M, V, [18] 7.0M, [17] 10.1[19] 5.9M, [18] 7.4V, [18] 5.7M, [17] NEST1 8.8BM 9.2M 8.8V 7.3BM 6.0M 7.3V 15.9BM 16.4M 15.9V 7.1BM 6.1M 7.1V N/A N/A PAGE 70 57 4.5 Discussion The equilibrium structures calculated w ith the vdW correction clearly show an improvement over pure DFT in the agreemen t with experiment. By pure DFT, the average error compared with the experimental data chosen for comparison is about 9%, and the vdW correction shows an improved error of approxi mately 2.3%. On the other hand, the results of Neumann and Perrin for 31 molecular crystals s howed an error of approximately 20% for pure DFT and about 1% with their vdW correction [2]. The difference could be the result of several changes that were made to the vdW correction for this work. It is expected that the cha nge of functional from PW91 to PBE would not contribute largely to the difference, as these functionals are commonly believed to provide very similar results. However, there were changes in this work for the energy cutoff, k point sampling density, and the tolerance for structural optimization that may be mostly responsible for the difference between the average results. The energy cutoff was 700 eV, yielding a larger basis than the 500 eV calculations of Neumann and Perrin. In addition, this work has a greater k point sampling density, which was approximately 0.05 1 near equilibrium conditions. The changes introduced for the energy cutoff and kpoint sampling in this work improve the accuracy of the DFT contribution. Also, the structural optimization in this work was set via a maximum force on any atom, but their work involved a tolerance of 0.010.02 kcal/mol in absolute energy. The different tolerance criteria could also contribute to the diffe rence because the ener gy tolerance of 0.010.02 kcal/mol (0.40.8 meV) for several system s was reached before the maximum force on any atom was less than 0.03 eV/. If re laxations were ceased at a total energy PAGE 71 58 convergence of 0.4 meV, the calculated unitcell volumes in most cases would be slightly greater. Further, the dampingfunction paramete rs of the vdW approach were not refit to adjust for the changes made in this work. It was expected that thes e parameters would not differ significantly with the changes made to the approach in the DFT component, but the contribution might be greater than expected. In either case, the vdW correction is a step towards improving the predictive po wer of DFT for these systems. For the hydrostaticcompression calcula tions, the vdW correction shows an improvement in the agreement with th e experimental EOS. While pure DFT overestimates the volume at a given pressure the vdW correction gi ves a slightly lower pressure at a given volume. One must consider that the experi mental data for most of the systems studied were taken at room temperat ure, but a comparison is made with DFT, a groundstate theory that gives results corresponding to 0 K. If thermal expansion on the 0 K results were accounted for, it is expected that the pressure at a given volume would increase. Hence, an accurate prediction at 0 K should predict a slightly smaller pressure at a given volume than an experiment at 300 K. Zerilli and Kuklja have provided an estimate of the pressure difference in the EM, FOX7, to account for the corresponding temperature difference of groundstate electronic structure calc ulations and experiment at 300 K [73], and the estimate at a density of 1900 kg/m3 was approximately 0.6 GPa. If a pressure increase on this order of magnitude were added to the vdWDFT results for the EOS, the agreement with most experime ntal studies would be much improved. The greatest contributing factor to the uncer tainty in the DFT calculations of the EOS is most likely the error in the pre ssure from the energy cutoff, which is approximately 0.4 GPa for the systems in this study. On the pressure scale of the EOS PAGE 72 59 plots, this uncertainty is roughly the same size as the sy mbols for the data points. Meanwhile, the error in the equilibriumst ructure calculations from the cutoff choice should also be considered. As shown in Ref [38], the equilibrium volume has a significant dependence on the basis, i.e. the en ergy cutoff, for several energetic molecular crystals. Convergence of unitcell volume with cutoff is notoriously slow for these systems, where the trend in approaching convergence is that volume increases with cutoff. Performing the calculations with a higher energy cutoff would decrease the error and most likely yield results for the vdW corr ection in better agreement with experiment, but this would incur a tremendous cost. Ov erall, the choice of the parameters was considered as a reasonable balance between accuracy and expense for the purposes of this study. It was expected that the influence of disp ersion is only significant at low pressure for molecular crystals. From th e results of the study, it is observed that the discrepancy between pure DFT and vdWDFT is greatest at low pressure, and the curves representing the EOS results for pure DFT and vdWDFT appear to approach one another with increasing pressure. However, the difference in pressure at a given volume is still significant at much higher pressures than expected. Although the damping function proposed by Neumann and Perrin is a breakthro ugh for equilibriumstructure calculation, further refinement of the damping function mi ght be necessary for calculations at high pressure. The bulk modulus predictions exhibit a few trends. The values calculated with the vdW correction are greater than those obtaine d with pure DFT. This result is expected since the vdW correction predicts a smaller volume for the unit cell of each system, and PAGE 73 60 the denser unit cell is most likely more resi stant to changes in volume. Further, the experiments of Gump and Peir is indicate that the bulk modul us decreases with increasing temperature [66]. If one were to account for the temperature difference between the simulations and experiment, it is possible that the vd W results would show better agreement with experimental values. Another trend exhibited is that the vdW correction appears to have little influence on the calcula ted value of the pressure derivative of the bulk modulus. An evaluation of the results obtained for the bulk modulus by comparison with reported experimental results is a difficult ta sk. The values from experiment seldom agree within an experimental uncertainty, and the discrepancy in the values measured (without uncertainty) is reflected in Table IV. This coul d be the result of seve ral factors, including the fitting forms used as well as a lack of lo wpressure data [1]. Fo r this reason, it is not clear whether the vdW correction improves or weakens the predictive ability of DFT in the prediction of these elastic properties. PAGE 74 61 [Manual chapter break] CHAPTER 5 UNIAXIALCOMPRESSION SIMULATIONS OF ENERGETIC MATERIALS Singlecrystal samples of explosives can respond quite differently to stimuli in comparison with polycrystalline samples. In the case of PETN, single crystals have shown greater initiation sensit ivity, and this behavior is also highly anisotropic. By investigating the factors that contribute to the anisotropic sensitivity of singlecrystal PETN, it might be possible to gain insight in to the poorly understood atomicscale events that lead to detonation. Following the discovery of unexpected shock response of PETN [91], the anisotropic sensitivity of PETN has been studied in several ex periments [23, 24, 9295] by Dick The initial works [23, 93] employed we dgetests that delive red plane shocks to samples of singlecrystal PETN in the <100> <001>, <110>, and <101> orientations (the notation < hkl > specifies a direction that is perpendicular to a cr ystallographic plane with similar Miller indices [93]). Reflected light from the free surface of the sample was used to determine the time after impact and distance into sample when and where, respectively, detonation took place. The resu lts show heavy directional dependence of sensitivity. The <110> and <001> directions were found to be sensitive. The <110> direction exhibits the greatest sensitivit y, where longitudinal st resses as low as approximately 4 GPa have caused detonation [93]. The <001> direction also has been described as sensitive because it has shown[ 23] a transition to detonation near 12 GPa. PAGE 75 62 On the other hand, the other two directions st udied in the experiments, <100> and <101>, are relatively insensitive. Th e initial works of Dick [23, 93] showed an intermediate velocity transition for the <101> direc tion, which indicated decomposition without detonation. The <100> directi on, however, showed no such transition [23, 93], and detonation was not observed for this direction until a later study by Yoo et al [25]. The shocktodetonation transition was achieved fo r <100> only near the detonation pressure of approximately 31 GPa [25]. Other energetic materials might also have anisotropic sensitivity properties, but there are no systems other than PETN for which the anisotropy has been sufficiently studied. Singlecrystal experiments of uni axial compression have been performed on HMX [96] and RDX [97], but the focus of the experiments was to obtain data on the elastic response to weak shocks. However, both works [96, 97] state that the anisotropic sensitivity of these materials will be investigated. Also, an isotropic sensitivity has been observed in a DAC experiment [98] on solid nitromethane, but the sensitive directions have not clearly been identified [99]. The mechanisms responsible for the anisot ropic sensitivity of PETN are not fully understood, though several works have suggested possibilities. The discoverer of the anisotropy, Dick, reasoned that directions with fewer available slip systems do not possess the ability to effectively relieve shea r stress, and these directions may be more sensitive than others [23, 24, 9294]. Comm only called the sterichindrance hypothesis, compression studies of PETN as well as other EMs have been designed based upon the available slip systems of the studied uni axialcompression directions [96, 97, 100]. Gruzdkov and Gupta [101] built upon the hypothe sis of steric hindrance by suggesting PAGE 76 63 that hindered shear in PETN causes confor mational changes in the PETN molecules which lead to local polarizati on of the crystalline lattice. The work assumes that these local changes can lead to ionic reactions, i.e. initia tion, under shock conditions. Another suggestion for the anisotropic sensitivity in P ETN was provided by Jindal and Dlott [102]. In their work, it is hypothesized that the anisotropic heating of the crystal as a result of uniaxial compression might be responsible for the observed directional dependence of sensitivity. Their results from an anharmonic potential for the molecular crystal naphthalene indicate that th e difference in temperature at pressures of approximately 10 GPa can lead to reaction rate constants that differ by orders of magnitude for different compression directions. Further, it has been proposed that the electronic band gap of EMs under compression might play a significa nt role in the initiation of detonation. In several works, Kuklja and coworkers [103107] have inve stigated the possible role of electronic excitations in the initiation of EMs. The effect of dislocations in the initiation of EMs is central to the work, where reduction of th e band gap by compression is a corequisite. Further, Reed et al. performed theoretical studies of nitromethane, concluding that dynamics outweigh static effects in the re duction of the bandgap [ 77] and a transient metallic layer is formed at the detonation front [108]. In this work, DFT is employed to inve stigate the response of EMs to uniaxial compression. Of particular intere st is the compression of PETN in the directions studied by in the experiments of Dick [23, 24]. Using the available experime ntal data regarding the anisotropic sensitivity of PETN, a comparison of th e physical properties for each PAGE 77 64 compression direction may be useful for id entifying physical quant ities that may be correlated with greater (or lesser) sensitivity. An aim of the study was to compare and contrast the shear stress during compression for several directions in each EM. The sterichindrance hypothesis emphasizes the role of shear in initiation, but other evidence also suggests its importance. For example, molecular dynamics simulations of uniaxial compression of RDX in the <100> direction by Cawkwell [100] show z ones of high shear as regions of energy localization. The melting of th e RDX occurs at 45 to the compression direction, where the shear stress is at a maximum. One may de duce that these regions of high shear stress have an associated localization of energy th at could, according to hot spot theory [109], serve to overcome energy barriers for hi ghly exothermic chemical reactions. In addition, the band gap under compression can be calculated from DFT. While the method described below captures only the static effects that contribute to the reduction of the band gap, it is important to investigate the anisotr opic nature of these static effects. This work investigates the band gap of each system under uniaxial compression to predict the an isotropic bandgap changes. Finally, simulations of EMs at the mesoscale require quantummechanical calculations of anisotropic constitutive re lationships. Menikoff and Sewell [110] have emphasized the need of anisotropic stressstr ain relationships for mesoscale simulations of HMX to accurately model th e concentration of stress between explosive grains, an effect that is believed to generate hot spots for detonation. In this work, these anisotropic constitutive relationships ar e predicted for several impor tant EMs for compression in PAGE 78 65 several lowindex crystallographic directions at pressures that ar e expected to be useful as input for simulation at greater length scales. Some of the results of this study have a ppeared in previous publications [79, 80, 111113] and appear below with permission where applicable. 5.1 Computational Details The uniaxial compressions were performed on the EMs in the lowindex crystallographic directions <100>, <010>, <001>, <110>, <101>, <011>, and <111>. It is important to note that the co mpression directions specified by are defined as perpendicular to crystallographic plan es referred to by Miller indices following the notation of Dick [24]. Clearl y, a direction specified by this notation is not equal (for general crystal symmetry) to the direction which in popular notation corresponds to the direction along where , and are the lattice vectors. To simulate uniaxial compression, the unit cell was first rotated such that the desired direction was aligned with the x axis. The compression was carried out in steps of 2% of the calculated equilibrium volume by first scaling the x components of the lattice vectors to reduce the volume of the unit cell by 2%, and then the atomic coordinates were relaxed. During the relaxations, the unitcell shape and atomic coordinates were held constant. Uniaxialcompression simulations were performed in each of the seven lowindex directions from 100% to 70% of th e calculated equilibrium volume. With the exception of the constraints during rela xation, the parameters for the VASP DFT PAGE 79 66 uniaxialcompression calculations are the same as those previously described on page 40 for hydrostatic compression. Note that th e uniaxialcompression calculations were performed without the vdW correct ion unless otherwise indicated. 5.2 Stressstrain relationships The constitutive relationships between stre ss and strain in the unit cells of the uniaxially compressed EMs of th is study are discussed below. 5.2.1 Principal stresses The stress tensor is used to describe the state of stre ss in a material. Owing to the fact that the elements of a tensor depend on the chosen coordinate axes, it is informative to extract invariant information on the stress state. This can be accomplished by diagonalizing the stress tensor to find its eigenvalues. The eigenvalues of the stress tensor are the principal stresses, and these quantitie s describe purely compressive stress on the corresponding orthogonal eigenvector s, i.e. the principal axes. The principal stresses have been defined in this work as , and where the subscripts are assigned in d ecreasing order. The greatest of the principal stresses corresponds to the longitudinal stress along the compression direction, while and correspond to the stresses along mutua lly perpendicula r directions. PAGE 80 67 The principal stresses , and of PETN under uniaxial compression are Figure 18: Principal stress as a function of volume ratio for seven lowindex uniaxial compressions. The hydrostatic pressure is shown for comparison. PAGE 81 68 shown in Figure 18, Figure 19, and Figur e 20, respectively. The principal stress clearly shows the anisotropic response of PETN to mechanical compression. There are directions with very similar behavior, and this res ponse is expected because of the crystal symmetry. In particular, the <100> and <010> compressions, as well as the <101> and <011> compressions, have almost identical principalstress values. The tetragonal symmetry of PETN indicates that the lattice constants a and b are equal, and one might expect these directions to exhibi t similar response to compression. Figure 19: Principal stress as a function of volume ratio. The pressure from the hydrostatic compressions is shown for comparison. PAGE 82 69 The greatest longitudinal stress for P ETN is observed at high compression along the <001> direction. The <110>, <101>, and <011> compressions show similar values at a volume ratio near 0.70. The smallest longit udinal stress is observed for the <100> and <010> directions. For the other pr incipal stresses, the values of are relatively isotropic with respect to the other stresses. However, the values of for the insensitive direction <100> and <010> are slightly smaller than the other directio ns of the study. Figure 20: Principal stress as a function of volume ratio for uniaxial compression in seven lowindex crystallographic directi ons of PETN. The pressure from the hydrostaticcompression simulations is shown for comparison. PAGE 83 70 The principal stresses , and of HMX are shown in Figure 21, Figure 22, and Figure 23, respectively. In contrast to the behavior shown by PETN, the longitudinal stresses from the uniaxialcompression simulations of HMX have very similar values at a volume ratio of 0.70. However, the <110> dir ection exhibits slightly greater values of than the other directions at high co mpression. Also, the < 011> direction shows relatively larger stress at lower compressi on, but nonmonotonic be havior is observed near a volume ratio of 0.80 that reduces the st ress to values that are similar to the other Figure 21: Principal stress (longitudinal stress) of HM X as a function of volume ratio for uniaxial compression in seven lo windex crystallographic directions. The pressure from the hydrostaticcompression calculations is sh own for comparison. PAGE 84 71 directions. With the exception of the increase in stress shown for the <100> and <001> directions near a volume ratio of 0.76, the values of for HMX are mostly isotropic. The values of principal stress are also relatively similar, where the <100> direction yields a slightly smaller stress th an the other compression directions. Figure 22: Principal stress of HMX as a function of volume ratio for uniaxial compression in seven lowindex crystallograp hic directions. The pressure from the hydrostaticcompression calculati ons is shown for comparison. PAGE 85 72 The principal stresses , and of RDX upon uniaxial compression are shown in Figure 24, Figure 25, and Figure 26, respectively. At high compression, the <011>, <100>, <001>, and <010> compressions show relatively greater stress. The values for RDX, like PETN and HMX, ar e relatively similar for the compression directions studied. There is a noticeable drop in the value of observed for the <110> direction near a volume ratio of 0.76. For the values of the <011> compression is Figure 23: Principal stress of HMX as a function of volume ratio for uniaxial compression in seven lowindex crystallograp hic directions. The pressure from the hydrostaticcompression calculati ons is shown for comparison. PAGE 86 73 significantly smaller than the values obtaine d for the other compression directions, which are relatively similar. Figure 24: Principal stress (longitudinal stress) of RDX as a function of volume ratio for uniaxial compression in seven lo windex crystallographic directions. The pressure from the hydrostaticcompression calculations is sh own for comparison. PAGE 87 74 Figure 25: Principal stress of RDX as a function of volume ratio for uniaxial compression in seven lowindex crystallograp hic directions. The pressure from the hydrostaticcompression calculati ons is shown for comparison. PAGE 88 75 The principal stresses , and of solid nitromethane from uniaxialFigure 26: Principal stress of RDX as a function of volume ratio for uniaxial compression in seven lowindex crystallograp hic directions. The pressure from the hydrostaticcompression calculati ons is shown for comparison. PAGE 89 76 compression simulations are shown in Figur e 27, Figure 28, and Figur e 29, respectively. The <001> direction stands out from the other directions as having the greatest value of the longitudinal stress As observed for the other systems, the values for the uniaxial compressions of nitrom ethane are relatively similar. On the other hand, the <010> compression has a significantly higher value of principal stress than the other compression directions. Figure 27: Principal stress (longitudinal stress) of nitromethane as a function of volume ratio for uniaxial compression in seve n lowindex crystall ographic directions. The pressure from the hydrostaticcompre ssion calculations is shown for comparison. PAGE 90 77 Figure 28: Principal stress of nitromethane as a function of volume ratio for uniaxial compression in seven lowindex cr ystallographic direct ions. The pressure from the hydrostaticcompression calculations is shown for comparison. PAGE 91 78 The principal stresses , and of NEST1 are shown in Figure 30, Figure 31, and Figure 32, respectively. Note that th e NEST1 uniaxialcompression calculations were performed with the vdW correction, but the other system s were treated with pure DFT. The values of for the uniaxial compressions of NEST1 are relatively similar for most of the compression range studied, but th e directions <010> and <110> display nonFigure 29: Principal stress of nitromethane as a function of volume ratio for uniaxial compression in seven lowindex cr ystallographic direct ions. The pressure from the hydrostaticcompression calculations is shown for comparison. PAGE 92 79 monotonic behavior near a volume ratio of 0.74. The <011> direction also shows a change in the trend of th e data that reduces the st ress near this volume ratio. Figure 30: Calculated maximum principal stress of NEST1 using vdWDFT for uniaxial compressions. The pressure from the hydrostaticcompression simulations is shown for comparison. PAGE 93 80 Figure 31: Calculated maximum principal stress of NEST1 using vdWDFT for uniaxial compressions. The pressure from the hydrostaticcompression simulations is shown for comparison. PAGE 94 81 5.2.2 Shear stresses The maximum shear stress is calculated from the principal stresses. The principal stresses are the eigenvalues of the stress tens or, and therefore they are the elements of a diagonalized stress tensor. By transforming th e stress tensor to find the axes where the nondiagonal elements, i.e. the shear stress es, are maximal, one discovers that the Figure 32: Calculated maximum principal stress of NEST1 using vdWDFT for uniaxial compressions. The pressure from the hydrostaticcompression simulations is shown for comparison. PAGE 95 82 maximum shear occurs near 45 to the princi pal axes, i.e. the eigenvectors corresponding to the principal stresses. To demonstrate the calculation of the sh ear stresses, the transformation of the stress tensor has the form (5.1) where are direction cosines and are stress tensor components in the transformed coordinates. Since the stress tensor with principalstress values is diagonal, Keeping the notation of the principal stresses, is denoted as principal stress As an example, one of the shear stresses is obtained by counterclockwise rotation of the x and y axes around the z axis by 45. By transformation, the maximum shear stress can be calculated as (5.2) The shear stress is equal in magnitude but opposite in sign for a clockwise rotation, hence only the magnitude of the shear stress is consid ered below. The quantities of interest from rotation of the principa l axes are therefore , and The greatest of these shear stresses is referred to as below, and the values of and are assigned in order of decreasing magnit ude. Note that the smallest shear stress maximum is typically much smaller than the ot her two. Owing to this fact, the values of for the uniaxialcompression simulations ar e expected to be of relatively minor PAGE 96 83 importance, but are shown below to give a more complete picture of the state of stress in the materials. The shear stresses , and for each uniaxial compression direction of PETN are shown in Figure 33, Figure 34, and Figur e 35, respectively. The <001> compression calculations show much greater shear stress than the other directions, which have similar values near a volume ratio of 0.70. The <101>, <011>, and <111> directions show nonmonotonic dependence of stress on strain below a volume ratio of 0.80. Figure 33: Maximum shear stress for uniaxial compressions of PETN. PAGE 97 84 Figure 34: Maximum shear stress for uniaxial compressions of PETN. PAGE 98 85 The ordering of the principal stresses by value masks the relation of the stress to crystal orientation, which is a very important consideration for PETN. To illustrate this point, the shear stresses and defined by (5.3) are shown in Figure 36 and Figur e 37, respectively. A few idea s should be considered for these stresses. First, the stress tensors obtai ned from the uniaxial compressions of PETN were approximately diag onal, so the values correspond to the principal stresses with orientational ordering. Several of the other sy stems did not yield ne arly diagonal stress Figure 35: Maximum shear stress for uniaxial compressions of PETN. PAGE 99 86 tensors, and the principalstress analysis was introduced to treat each system in an equal Figure 36: Shear stress from uniaxialcompression simulations for directions of known sensitivity in PETN. PAGE 100 87 manner. Next, note that the compressi on direction in all cases is the x direction, and therefore the shear stresses and reflect the maximum shear stress at 45 to the compression direction ( is considered to be of minor importance). In Figure 36 and Figure 37, these maximum shear stresses ar e shown for the directions of known sensitivity in PETN. Recall that the <110> a nd <001> directions are sensitive, while the <100> and <101> directions are insens itive to shock compression. Although the values do not appear to largely differ between the high and low sensitivity directions, Figure 37: Shear stress from uniaxialcompression simulations for directions of known sensitivity in PETN. PAGE 101 88 there is a clear distinction in the values of between the sensitive and insensitive directions. At a volume ratio of 0.70, the insens itive directions exhibit much lower values of shear stress than the sensitive directions. Furthe r, both of the insensitive directions have nonmonotonic dependence of shear stress on volume ratio. These observations indicate that there might be a correlation between shock sensitivity and maximum shear stresses in EMs. The maximum shear stresses , and from the uniaxialcompression simulations of HMX are shown in Figur e 38, Figure 39, and Figure 40, respectively. Figure 38: Maximum shear stress from uniaxialcompression simulations of HMX. PAGE 102 89 Compression in the <110> directio n yields large values of both and According to the possible correlation of sensitivity with gr eater shear stress observed in PETN, an investigation into th e relative sensitivity of the < 110> direction in HMX would be enlightening. Also, all of the directions show nonmonotonic dependence in one or both of the larger two shear stresses except for the <010> and <001> directions. In PETN, nonmonotonicity of the shear stresses as a f unction of volume ratio was observed for the insensitive directions, and this behavior might also re flect insensitivity in HMX. Figure 39: Maximum shear stress from uniaxialcompression simulations of HMX. PAGE 103 90 Figure 40: Maximum shear stress from uniaxialcompression simulations of HMX. PAGE 104 91 The maximum shear stresses , and from the uniaxialcompression simulations of RDX are shown inFigure 41, Figure 42, and Figure 43, respectively. The directions with greater shear stress, and possibly greater sensitivity, at high compression are <100>, <010>, <001>, and <011> Als o, the <101> and <111> directions show much smaller shear stresses upon compre ssion, and therefore might be relatively insensitive directions. Figure 41: Maximum shear stress from uniaxialcompression simulations of RDX. PAGE 105 92 Figure 42: Maximum shear stress from uniaxialcompression simulations of RDX. PAGE 106 93 Figure 43: Maximum shear stress from uniaxialcompression simulations of RDX. PAGE 107 94 The maximum shear stresses , and from the uniaxialcompression simulations of nitromethane are show n in Figure 44, Figure 45, and Figure 46, respectively. The directions that stand out fr om the others are the <001> direction, which shows much greater shear stresses and and the <010> direction, which yields smaller shear stresses than the other compressi on directions. The othe r directions reveal relatively similar behavior. Figure 44: Maximum shear stress from uniaxialcompression simulations of nitromethane. PAGE 108 95 Figure 45: Maximum shear stress from uniaxialcompression simulations of nitromethane. PAGE 109 96 Figure 46: Maximum shear stress from uniaxialcompression simulations of nitromethane PAGE 110 97 The maximum shear stresses , and from the uniaxialcompression simulations of NEST1 are shown in Fi gure 47, Figure 48, and Figure 49, respectively. All directions show nonmonotonic behavior, but the <100> and <001> directions show greater shear stress values for and at high compression. Figure 47: Maximum shear stress for uniaxial compressions of NEST1 using vdWDFT. PAGE 111 98 Figure 48: Maximum shear stress for uniaxial compressions of NEST1 using vdWDFT. PAGE 112 99 5.3 Band Gaps As discussed above, the reduction of the electronic band gap in EMs during compression has been identified as a possible contributing factor to initiation. From the compression simulations, the band gap was approximately calculated to determine if the gap was reduced significantly and to estimat e the anisotropy in the bandgap reduction caused by compression. Figure 49: Maximum shear stress for uniaxial compressions of NEST1 using vdWDFT. PAGE 113 100 For the band gap cal culations, only the k points used in the VASP calculations were included. It is expected that the few points used fo r sampling provide a reasonable estimate of the band gap. In ot her words, it is assumed that there are not la rge variations in the band structure for these materials. The band gaps as functions of volum e ratio for the uniaxialcompression simulations of PETN are shown in Figure 50. From a comparison with the band gap from the hydrostaticcompression simulations, the uni axial compressions lower the calculated Figure 50: Band gap of PETN from uniaxi alcompression simulations. The band gap from the hydrostaticcompression simu lations is shown for comparison. PAGE 114 101 band gap of PETN more than hydrostatic comp ression. However, it is clear that none of the band gaps approach a value typical of a conductor. In additi on, the calculated band gap does not appear to correlate with sens itivity. The most sensitive direction, <110>, shows the least change in band gap, and the least sensitive direc tion, <100>, has only a slightly greater reduction. Fu rther, the sensitive directi on <001> has a greater reduction in band gap than <100>, but the insensitive di rection <101> has the greatest change in band gap of the directions in the study. PAGE 115 102 The band gaps as functions of volum e ratio for the uniaxialcompression simulations of HMX are shown in Figure 51. The <011> direction displays a much greater reduction in the band gap than the ot her directions. On the other hand, the <001> and <101> directions show less of a re duction in the band gap than hydrostatic compression. While the calculations indicate a significant anisotropy in the band gap for HMX, it does not approach a metallic state. Figure 51: Band gap of HMX from uniaxia lcompression simulations. The band gap from the hydrostaticcompression simu lations is shown for comparison. PAGE 116 103 The band gaps as functions of volum e ratio for the uniaxialcompression simulations of RDX are shown in Figure 52. The <100> direction exhibits a greater reduction than the other directions. Als o, the <011>, <010>, and <001> compressions display a greater gap than hydrostatic comp ression at intervals within the compression range shown. Note that the band gap at equilib rium conditions is similar to that of HMX, and the minimum gap at high compressi on is also approximately equal. Figure 52: Band gap of RDX from uniaxi alcompression simulations. The band gap from the hydrostaticcompression simu lations is shown for comparison. PAGE 117 104 The band gaps as functions of volum e ratio for the uniaxialcompression simulations of nitromethane are shown in Figure 53. Note the increase in band gap from equilibrium to the first compression step for so me directions. This behavior is assumed to be caused by the error in the approxima tion used to calculate the band gap. The equilibrium calculations we re performed with the k points used for the hydrostatic compressions, but each uniaxial compression had different points used for sampling. Hence, the increase in band gap is mo st likely a result of the change in k points, and the energy difference between the band gap from the hydrostatic calculations at a volume Figure 53: Band gap of nitromethane fr om uniaxialcompression simulations. The band gap from the hydrostaticcompression simulations is shown for comparison. PAGE 118 105 ratio of 0.98 and the <010>, <001>, and <011> results is a reasonable estimate of the minimum error in the approximation used to calculate the band gap. The band gaps as functions of volum e ratio for the uniaxialcompression simulations of NEST1 are shown in Figur e 54. At high compression, the <110> direction exhibits a relatively greater reduction in the band gap as compared to the other directions. The <101> and<111> directions show an unexp ected increase in the band gap near a volume ratio of 0.80. For all directions, the reduction of the band gap is below 1 eV. Figure 54: Band gap of NEST1 from uni axialcompression simulations. The band gap from the hydrostaticcompression simulations is shown for comparison. PAGE 119 106 5.4 Discussion The principal and shear stresses, as well as the band gap, as a function of volume ratio were examined for uniaxial compressions of PETN, HMX, RDX, nitromethane, and NEST1. Because the anisotropic shock se nsitivity of PETN has been studied by experiment, the response of PETN was examin ed for correlations with sensitivity. It was observed that the insensitiv e <100> and <101> directions in PETN exhibited lower maximum shear stresses, as well as nonmonotonic dependence on strain, and the sensitive <110> and <001> directions clea rly showed greater maximum shear stress values. Although further work is needed to validate a correlation, the results present convincing evidence that a correlation is possi ble. If, in fact, there is a correlation and it can be extended to other EMs, it would be valu able to identify direc tions in other systems that display similar behavior to the sensitiv e and insensitive direc tions in PETN. Hence, the anisotropic shear stresses of HMX, RDX, nitromethane, and NEST1 were calculated, and directions with relative ly greater (and lesser) shear stress were identified. The impending results of experiments by other groups [96, 97] for these materials will provide a means to accept or reject the possible correlation. The band gap of the materials was also examined for a correlation with sensitivity in PETN, but it was found to be unlikely fr om the results. Meanwhile, the EMs of the study exhibit anisotropic reduction of the band gap with compression. Overall, the calculated band gaps reduced by roughly 1 eV, clearly not enough to cause metallization. As discussed above, the simulations do not pr ovide a complete picture for reductions in PAGE 120 107 the band gap because of the static nature of the calculations [77] and because the calculations involve defe ctfree unit cells [104]. The uniaxialcompression calculations performed on NEST1 differ from the other systems because they were performed with the vdW correction. Calculations have been performed to compare th e results between pure DFT a nd the vdW correction for the principal and shear stresses. A preliminary insp ection of the results indicate that the vdW correction gives higher values for the prin cipal stresses (and, consequently, the shear stresses), but the relative differences in st resses between the comp ression directions are similar. Hence, it is expected that directions showing greater princi pal and/or shear stress from pure DFT calculations would also exhibi t greater relative values with the vdW correction. Meanwhile, the pure DFT calculations for all systems were not repeated with the vdW correction for a thorough comparis on owing to the large expense of the calculations and because a si gnificant difference in the relative behavior was not expected. In the stresses and band gaps for some uniaxial compressions, abrupt changes are observed, such as the onset of nonmonotonic be havior and sudden large increases. While an analysis of the changes in molecular structure that accompany these interesting changes would be valuable, this investiga tion has not yet been thoroughly performed. However, it is expected that those in co llaboration with this study will examine how changes to molecular structure influence the electronic structure and stress in these EMs. PAGE 121 108 [Manual chapter break] CHAPTER 6 CONCLUSION In this work, density functional theory was used to study the energetic materials (EMs) PETNI, HMX, RDX, solid nitromethane, and a recently discovered EM referred to as NEST1. The empirical van der Waals correction of Neumann and Perrin [2] was used along with the DFT code VASP in an effort to improve the description of the molecular crystals of the study. The empi rical method was used in the prediction of equilibrium unitcell structures and h ydrostaticcompression simulations, and a comparison was made between the corrected and uncorrected results. For the prediction of equilibrium unitcell volumes, the pure DFT results provided an error of about 9% when averaged over each of the five EMs of the study, whereas the error was reduced to approximately 2% w ith the vdW correction. Although this is a significant improvement, the average percent er ror from the work of Neumann and Perrin is greater for the pure DFT results and smaller for the corrected results. It was concluded that a large contribution to the discrepancy was likely caused by the changes in this study to the DFT portion of the empirical method. Th e changes were introduc ed with the intent to improve the quality of the DFT portion. A worthwhile study would be to repeat the fitting of the vdW potential parameters with altered VASP settings to investigate whether a similar agreement of 1% is obtained [2]. PAGE 122 109 The isothermal equations of state (EOSs) at 0 K for each EM were calculated from the hydrostaticcompression simulations a nd compared with available experimental data. For PETN and RDX, there is a cl ear improvement in the agreement with experiment. For HMX and nitromethane, there ar e data from two or more experiments for comparison, and the agreement is similar to, if not better than, the agreement yielded by pure DFT. Meanwhile, the vdW correction predic ts a smaller volume at a given pressure, which is physically consistent with th e difference between the simulation and experimental temperatures. Note that expe rimental data on the hydrostatic compression of NEST1 are not currently available for comparison. The calculated lattice constants a b and c from the hydrostaticcompression simulations were also compared with experi ment. With the exception of a few cases, the vdW correction increases the agreement with experiment. The correction therefore shows an improvement in the agreement on average for the EMs studied. The bulk modulus and its pressure derivati ve were calculated from the hydrostatic compression data and compared with availabl e experimental data. The BirchMurnaghan, Murnaghan, and Vinet EOSs were used to f it the simulation data. Owing to the lack of agreement in experimental results, it was difficult to make a sound conclusion regarding the improvement of the vdW correction in th e prediction of the bulk moduli. However, trends were observed for the re sults; the pure DFT approach predicted a smaller value for the bulk modulus than experiment in ma ny cases, whereas the vdWDFT approach yielded greater values than experiment on average. PAGE 123 110 For the comparison, the temperature differe nce must be considered between the simulations and the experimental data. The experimental data were taken at room temperature, but the groundstate results repr esent the conditions at absolute zero. Thus, it is expected that the theoretical results should predict a smaller volume (and latticevector lengths, correspondingl y) at a given pressure than experiment. The vdWcorrection results are consistent with this expectation, whereas the pure DFT results are not. Meanwhile, it remains to be shown th at the difference in volume between the experimental data and the vdWDFT results at a given pressu re is consistent with the temperature difference. A suggestion for future research is the addition of methods to account for temperature to the simulati ons, which would allow a more suitable comparison with experiment and enable furthe r exploration into the equation of state for EMs. Uniaxialcompression simulations were also performed in seven lowindex crystallographic directions fo r the EMs of this study. The cal culations were designed to investigate the anisotropic response of th e materials to compression, with particular interest in the shocksensitive directions of PETN. The constitutive relationships between st ress and strain were calculated for the EMs under uniaxial compression. Specifically, the principal stresses and the maximum shear stresses were calculated for unia xial compression up to a volume ratio of V/V0=0.70. The predicted behavior of stress for e ach material is clearly anisotropic. Unlike the other EMs of the study, the anis otropic sensitivity of PETN to shock has been studied thoroughly by experiment [ 23, 24]. From a comparison of the maximum PAGE 124 111 shear stresses at about 45 to the compressi on direction, greater sh ear stresses were found for the sensitive <110> and <001> directions than for the insensitive <100> and <101> directions. Further, both in sensitive directions exhibi ted nonmonotonic dependence of shear stress on strain, and this behavior was not shown by the sensitive directions. It was concluded from these observations that th e maximum shear stresses might have a correlation with shock sensitivity. For this r eason, the directions co rresponding to greater and lesser shear stress es, as well as nonmonotonic depe ndence of shear stress on strain, were identified in HMX, RDX, nitromethane and NEST1. The results of experiments on the sensitivity of these materials is highly anticipated to determine if the sensitive (insensitive) directions from e xperiment correspond to direc tions with greater (lesser) shear stress from the simulations. The electronic band gap was al so calculated for each of the uniaxially compressed EMs. These results also indicate anisotropy, but the calculated band ga ps do not appear to be correlated with the anisotropic sensitivity in PETN. Further, the band gap reduces by approximately 1 eV at most for each EM at a volume ratio up to 0.70, which is not nearly enough to render the systems metallic. Owing to the expectation that dynamic effects and crystalline defects might be responsible for greater bandgap reduction [77, 103], it was anticipated that the band gaps calculated in this study would not decrease significantly. The stresses and band gaps as functions of volume ratio exhibited abrupt changes with uniaxial compression for some directi ons. The relation of th ese changes to the underlying molecular structure is not analyzed in this work, but this investigation will be the subject of future research. PAGE 125 112 REFERENCES [1] R. Menikoff, and T. D. Sewell, "Fitting forms for isothermal data", High Pressure Research 21 121 (2001). [2] M. A. Neumann, and M. A. Perrin, "Ene rgy ranking of molecular crystals using density functional theo ry calculations and an empirical van der Waals correction", Journal of Physical Chemistry B 109 15531 (2005). [3] H. H. Cady, and A. C. Larson, "Pentaer ythritol Tetranitrate Ii Its CrystalStructure and Transformation to Petn I an Algorithm for Refinement of CrystalStructures with Poor Data", Acta Crys tallographica Section BStructural Science 31 1864 (1975). [4] T. R. Gibbs, and A. Popolato, LASL Explosive Property Data (University of California Press, Berkeley, 1980). [5] O. Tschauner et al. "Structural transition of PETN I to ferroelastic orthorhombic phase PETNIII at elevated pressure s", Journal of Chemical Physics 127 (2007). [6] A. D. Booth, and F. J. Llewellyn, "T he Crystal Structure of Pentaerythritol Tetranitrate", Journal of the Chemical Society, 837 (1947). [7] J. Akhavan, The Chemistry of Explosives (The Royal Society of Chemistry, Cambridge, 2004). PAGE 126 113 [8] C. S. Choi, and H. P. Boutin, "A Study of Crystal Structure of BetaCyclotetramethylene Tetranitramine by Neut ron Diffraction", Acta Crystallographica Section BStructural Crystall ography and Crystal Chemistry B 26 1235 (1970). [9] J. A. Ciezak et al. "Highpressure vibr ational spectroscopy of energetic materials: Hexahydro1,3,5trinitro1,3,5triazine", Journal of Physical Chemistry A 111 59 (2007). [10] A. J. Davidson et al. "Explosives under pressure the crystal structure of gammaRDX as determined by highpressure Xra y and neutron diffrac tion", Crystengcomm 10 162 (2008). [11] Z. A. Dreger, and Y. M. Gupta, "High pressure Raman spectroscopy of single crystals of hexahydro1,3,5trin itro1,3,5triazine (RDX)", Jour nal of Physical Chemistry B 111 3893 (2007). [12] B. Olinger, P. M. Halleck, and H. H. Cady, "Isothermal Linear and Volume Compression of Pentaerythritol Tetranitrate (Petn) to 10 Gpa (100 Kbar and Calculated Shock Compression", Journal of Chemical Physics 62 4480 (1975). [13] B. Olinger, B. Roof, and H. Cady, in Commissariat a l'Energie Atomique Saclay, France, 1978), pp. 3. [14] C.S. Yoo et al. "Equations of State of Unr eacted High Explosives at High Pressures", in 11th International De tonation Synmposium (unpublished, Snowmass Village, CO, 1998). [15] C. S. Choi, and E. Prince, "CrystalStr ucture of CyclotrimethyleneTrinitramine", Acta Crystallographica Section BStructur al Crystallography and Crystal Chemistry B 28 2857 (1972). PAGE 127 114 [16] S. F. Trevino, E. Prince, and C. R. H ubbard, "Refinement of th e Structure of Solid Nitromethane", Journal of Chemical Physics 73 2996 (1980). [17] D. T. Cromer, R. R. Ryan, and D. Sc hiferl, "The Structure of Nitromethane at Pressures of 0.3 to 6.0 Gpa", J ournal of Physical Chemistry 89 2315 (1985). [18] M. Citroni et al. "Crystal structure of nitromethane up to the reaction threshold pressure", Journal of Physical Chemistry B 112 1095 (2008). [19] F. L. Yarger, and B. Olinger, "Compr ession of Solid Nitromethane to 15 Gpa at 298 K", Journal of Chemical Physics 85 1534 (1986). [20] S. Courtecuisse et al. "PhaseTransitions and ChemicalTransformations of Nitromethane up to 350DegreesC and 35Gpa", Journal of Chemical Physics 102 968 (1995). [21] R. Ouillon et al. "Raman and infrared investigati ons at room temperature of the internal modes behaviour in solid nitrometha neh(3) and d(3) up to 45 GPa", Journal of Raman Spectroscopy 39 354 (2008). [22] D. E. Chavez et al. "Synthesis of an Energeti c Nitrate Ester", Angewandte ChemieInternational Edition 47 8307 (2008). [23] J. J. Dick, "Effect of Crystal Or ientation on Shock Initiation Sensitivity of Pentaerythritol Tetranitrate Expl osive", Applied Physics Letters 44 859 (1984). [24] J. J. Dick, "Anomalous shock initiation of detonation in pentaer ythritol tetranitrate crystal", Journal of Applied Physics 81 601 (1997). [25] C. S. Yoo et al. "Anisotropic shock sensitivity and detonation temperature of pentaerythritol tetranitr ate single crystal", Jour nal of Applied Physics 88 70 (2000). PAGE 128 115 [26] M. Born, and R. Oppenheimer, "Qua ntum theory of molecules", Annalen Der Physik 84 0457 (1927). [27] R. M. Martin, Electronic Structure: Basic Theory and Practical Methods (Cambridge University Press, Cambridge, 2004). [28] R. Car, "Introduction to densityf unctional theory and abinitio molecular dynamics", Quantitative StructureActivity Relationships 21 97 (2002). [29] P. Hohenberg, and W. Kohn, "Inhomogene ous Electron Gas", Physical Review B 136 B864 (1964). [30] W. Kohn, and L. J. Sham, "SelfCons istent Equations Including Exchange and Correlation Effects", Physical Review 140 1133 (1965). [31] I. G. Kaplan, Intermolecular Interactions: Phys ical Picture, Computational Methods and Model Potentials (John Wiley & Sons, Ltd., 2006). [32] M. D. Segall et al. "Firstprinciples simulation: ideas, illustrations and the CASTEP code", Journal of PhysicsCondensed Matter 14 2717 (2002). [33] P. Gianozzi, in Density Functional Theory for Electronic Structure Calculations 2005). [34] J. Moreno, and J. M. Soler, "Optimal Meshes for Integrals in RealSpace and ReciprocalSpace Unit Cells", Physical Review B 45 13891 (1992). [35] H. J. Monkhorst, and J. D. Pack, "Speci al Points for BrillouinZone Integrations", Physical Review B 13 5188 (1976). [36] R. P. Feynman, "Forces in molecules", Physical Review 56 340 (1939). [37] G. Kresse, and J. Furthmuller, "Effici ent iterative schemes for ab initio totalenergy calculations using a planewa ve basis set", Physical Review B 54 11169 (1996). PAGE 129 116 [38] E. F. C. Byrd, G. E. Scuseria, and C. F. Chabalowski, "An ab initio study of solid nitromethane, HMX, RDX, and CL20: Succes ses and failures of DFT", Journal of Physical Chemistry B 108 13100 (2004). [39] Y. Andersson et al. "Densityfunctional account of van der Waals forces between parallel surfaces", Soli d State Communications 106 235 (1998). [40] E. Hult et al. "Density functional for van der Waals forces at surfaces", Physical Review Letters 77 2029 (1996). [41] E. Hult et al. "Unified treatment of asymptotic van der Waals forces", Physical Review B 59 4708 (1999). [42] W. Kohn, Y. Meir, and D. E. Maka rov, "vanderWaals energies in density functional theory", Phys ical Review Letters 80 4153 (1998). [43] K. Rapcewicz, and N. W. Ashcroft, "F luctuation Attraction in Condensed Matter a Nonlocal FunctionalAppro ach", Physical Review B 44 4032 (1991). [44] J. F. Dobson, and B. P. Dinte, "Con straint satisfaction in local and gradient susceptibility approximations: Application to a van der Waals density functional", Physical Review Letters 76 1780 (1996). [45] B. I. Lundqvist et al. "DensityFunctional Theory Including VanDerWaals Forces", International Journal of Quantum Chemistry 56 247 (1995). [46] M. Dion et al. "Van der Waals density functi onal for general geometries", Physical Review Letters 92 (2004). [47] M. Dion et al. "Van der waals density functiona l for general geometries (vol 92, art no 246401, 2004)", Physical Review Letters 95 (2005). PAGE 130 117 [48] D. C. Langreth et al. "Van der Waals density functional theory with applications", International Journal of Quantum Chemistry 101 599 (2005). [49] D. C. Langreth et al. "A density functional for sparse matter", Journal of PhysicsCondensed Matter 21 (2009). [50] H. Rydberg et al. "Van der Waals density functi onal for layered structures", Physical Review Letters 91 (2003). [51] H. Rydberg et al. "Tractable nonlocal correlati on density functionals for flat surfaces and slabs", Physical Review B 62 6997 (2000). [52] G. Kresse, and J. Furthmuller, "Effici ency of abinitio total energy calculations for metals and semiconductors using a planewave basis set", Computational Materials Science 6 15 (1996). [53] J. P. Perdew et al. "Atoms, Molecules, Solids, and Surfaces Applications of the Generalized Gradient Approximation for Excha nge and Correlation", Physical Review B 46 6671 (1992). [54] J. P. Perdew et al. "Atoms, Molecules, Solids, and Surfaces Applications of the Generalized Gradient Approximation for Ex change and Correlation (Vol 46, Pg 6671, 1992)", Physical Review B 48 4978 (1993). [55] J. P. Perdew, K. Burke, and M. Ernz erhof, "Generalized gr adient approximation made simple (vol 77, pg 3865, 1996)", Physical Review Letters 78 1396 (1997). [56] W. T. M. Mooij et al. "Transferable ab initio intermolecular potentials. 1. Derivation from methanol dimer and trimer cal culations", Journal of Physical Chemistry A 103 9872 (1999). PAGE 131 118 [57] Q. Wu, and W. T. Yang, "Empirical correction to density functional theory for van der Waals interactions", Journal of Chemical Physics 116 515 (2002). [58] T. A. Halgren, "Representation of Va nderwaals (Vdw) Interactions in Molecular Mechanics ForceFields Potential Form, Combination Rules, and Vdw Parameters", Journal of the American Chemical Society 114 7827 (1992). [59] A. Kumar, and W. J. Meath, "Pseudo Spectral Dipole OscillatorStrengths and DipoleDipole and TripleDipole Dispersion En ergy Coefficients for Hf, Hcl, Hbr, He, Ne, Ar, Kr and Xe", Molecular Physics 54 823 (1985). [60] A. Kumar, and W. J. Meath, "Di pole Oscillator Strength Properties and Dispersion Energies for Acetylene and Benzene", Molecular Physics 75 311 (1992). [61] G. R. Burton et al. "Valence Shell Absolute Ph otoabsorption OscillatorStrengths, Constrained Dipole Oscillator Strength Distri butions, and Dipole Properties for Ch3nh2, (Ch3)(2)Nh, and (C h3)(3)N", Canadian Journal of ChemistryRevue Canadienne De Chimie 72 529 (1994). [62] A. Kumar, and W. J. Meath, "Isotropic dipole properties for acetone, acetaldehyde and formaldehyde", Molecular Physics 90 389 (1997). [63] B. L. Jhanwar, and W. J. Meath, "Di pole Oscillator Strength Distributions, Sums, and Dispersion Energy Coefficients for Co and Co2", Chemical Physics 67 185 (1982). [64] B. L. Jhanwar, and W. J. Meath, "PseudoSpectral Dipole Oscillator Strength Distributions for the Normal Alkanes th rough Octane and the Evaluation of Some Related DipoleDipole and TripleDipole Di spersion Interaction Energy Coefficients", Molecular Physics 41 1061 (1980). PAGE 132 119 [65] C. S. Yoo, and H. Cynn, "Equation of state, phase transition, decomposition of betaHMX (octahydro1,3,5,7tetr anitro1,3,5,7tetrazo cine) at high pressures", Journal of Chemical Physics 111 10229 (1999). [66] J. C. Gump, and S. M. Peiris, "Isoth ermal equations of state of beta octahydro1,3,5,7tetranitro1,3,5,7tetrazocine at high temperatures", Jo urnal of Applied Physics 97 (2005). [67] C. K. Gan, T. D. Sewell, and M. Challacomb, "Allelectron densityfunctional studies of hydrostatic comp ression of pentaerythritol te tranitrate C(CH2ONO2)(4)", Physical Review B 69 (2004). [68] E. F. C. Byrd, and B. M. Rice, "Ab initio study of compressed 1,3,5,7tetranitro1,3,5,7tetraazacyclooctane (HMX), cyclot rimethylenetrinitramine (RDX), 2,4,6,8,10,12hexanitrohexaazaisowurzitane (CL20), 2,4,6trinitro1,3,5benzenetriamine (TATB), and pentaerythritol tetranitrate (PETN)", Journal of Physical Chemistry C 111 2787 (2007). [69] H. V. Brand, "Ab initio allelectr on periodic HartreeFo ck study of hydrostatic compression of pentaerythrito l tetranitrate", Journal of Physical Chemistry B 109 13668 (2005). [70] T. D. Sewell, "Monte Carlo simula tion of the hydrostatic compression of RDX and betaHMX." Abstracts of Papers of the American Chemical Society 215 U560 (1998). [71] D. C. Sorescu, B. M. Rice, and D. L. Thompson, "Theoretical studies of the hydrostatic compression of RDX, HMX, HNIW, a nd PETN crystals", Jo urnal of Physical Chemistry B 103 6783 (1999). PAGE 133 120 [72] T. D. Sewell et al. "A molecular dynamics simula tion study of elastic properties of HMX", Journal of Chemical Physics 119 7417 (2003). [73] F. J. Zerilli, and M. M. Kuklja, "Fir st principles calculation of the mechanical compression of two organic molecular crysta ls", Journal of Physical Chemistry A 110 5173 (2006). [74] D. Lian et al. "Highpressure behaviour of betaHMX crystal studied by DFTLDA", Chinese Physics Letters 25 (2008). [75] F. J. Zerilli, J. P. Hooper, and M. M. Kuklja, "Ab initio studies of crystalline nitromethane under high pressure", Journal of Chemical Physics 126 (2007). [76] H. Liu et al. "Structural and vibra tional properties of so lid nitromethane under high pressure by density functional theo ry", Journal of Chemical Physics 124 (2006). [77] E. J. Reed, J. D. Joannopoulos, and L. E. Fried, "Electronic excitations in shocked nitromethane", Physical Review B 62 16500 (2000). [78] M. W. Conroy et al. "Application of van der Waals density functional theory to study physical properties of energetic mate rials", in 16th APS Topical Conference on Shock Compression of Condensed Matter edited by M. L. Elert et al. Nashville, TN, 2009). [79] M. W. Conroy et al. "Firstprinciples studies of hydrostatic and uniaxial compression of a new energetic material an energetic nitrate ester", in 16th APS Topical Conference on Shock Compression of C ondensed Matter, edited by M. L. Elert et al. (AIP, Nashville, TN, 2009). PAGE 134 121 [80] M. W. Conroy et al. "Density Functional Theory Calculations of Solid Nitromethane under Hydrostatic and Uniaxi al Compressions with Empirical van der Waals Correction", Journal of Physical Chemistry A 113 3610 (2009). [81] J. P. Perdew, K. Burke, and Y. Wa ng, "Generalized gradient approximation for the exchangecorrelation hole of a manyelectron system", Physical Review B 54 16533 (1996). [82] P. E. Blochl, "Projector Augmen tedWave Method", Physical Review B 50 17953 (1994). [83] G. Kresse, and D. Joubert, "From ul trasoft pseudopotentials to the projector augmentedwave method" Physical Review B 59 1758 (1999). [84] G. Kresse, and J. Hafner, "NormCo nserving and Ultrasof t Pseudopotentials for FirstRow and TransitionElements", Journal of PhysicsCondensed Matter 6 8245 (1994). [85] D. Vanderbilt, "Soft SelfConsiste nt Pseudopotentials in a Generalized Eigenvalue Formalism", Physical Review B 41 7892 (1990). [86] P. R. Eiland, and R. Pepinsky, Z. Kristallogr. 106 (1955). [87] J. M. Winey, and Y. M. Gupta, "Sec ondorder elastic constant s for pentaerythritol tetranitrate single crystals", Journal of Applied Physics 90 1669 (2001). [88] S. Haussuhl, "Elastic and thermoelasti c properties of select ed organic crystals: acenaphthene, transazobenzene, benzophenone, tolane, transstilbene, dibenzyl, diphenyl sulfone, 2,2 'biphenol, urea, melamine, hexogen, succinimide, pentaerythritol, urotropine, malonic acid, dimethyl malonic acid, maleic acid, hippuric acid, aluminium PAGE 135 122 acetylacetonate, iron acetylacetonate, and tetraphenyl silicon", Zeitschrift Fur Kristallographie 216 339 (2001). [89] R. B. Schwarz et al. "Resonant ultrasound spectroscopy measurement of the elastic constants of cyclotrimethylene tr initramine", Journal of Applied Physics 98 (2005). [90] J. J. Haycraft, L. L. Stevens, and C. J. Eckhardt, "The elastic constants and related properties of the energetic material cyclot rimethylene trinitramine (RDX) determined by Brillouin scattering", Journal of Chemical Physics 124 (2006). [91] P. M. Halleck, and J. Wackerle, "Dyna mic ElasticPlastic Properties of SingleCrystal Petn", Bulletin of the American Physical Society 20 20 (1975). [92] J. J. Dick, "Supercritical Shear in Sh ocked Pentaerythritol Tetranitrate", Applied Physics Letters 60 2494 (1992). [93] J. J. Dick et al. "Shock Response of Pentaerythritol Tetranitrate SingleCrystals", Journal of Applied Physics 70 3572 (1991). [94] J. J. Dick, and J. P. Ritchie, "The Crystal Orientation Dependence of the Elastic Precursor Shock Strength in the Petn Molecu lar Explosive and the M odeling of the Steric Hindrance to Shear by Molecular Mech anics", Journal De Physique Iv 4 393 (1994). [95] J. J. Dick, and J. P. Ritchie, "Mol ecular Mechanics Mode ling of Shear and the Crystal Orientation Dependence of the Elastic Precursor Shock Strength in Pentaerythritol Tetranitrate", Journal of Applied Physics 76 2726 (1994). [96] J. J. Dick et al. "Elasticplastic wave prof iles in cyclotetramethylene tetranitramine crystals", Journal of Applied Physics 96 374 (2004). PAGE 136 123 [97] D. E. Hooks, K. J. Ramos, and A. R. Martinez, "Elastic plastic shock wave profiles in oriented single cr ystals of cyclotrimethylene trinitramine (RDX) at 2.25 GPa", Journal of Applied Physics 100 (2006). [98] G. J. Piermarini, S. Block, and P. J. Miller, "Effects of Pressure on the ThermalDecomposition Kinetics and ChemicalReactivity of Nitromethane", Journal of Physical Chemistry 93 457 (1989). [99] J. J. Dick, "OrientationDependent Expl osion Sensitivity of Solid Nitromethane", Journal of Physical Chemistry 97 6193 (1993). [100] M. J. Cawkwell et al. "Shockinduced shear bands in an energetic molecular crystal: Application of shockfront absorbi ng boundary conditions to molecular dynamics simulations", Physical Review B 78 (2008). [101] Y. A. Gruzdkov, and Y. M. Gupta, "Shock wave initiation of pentaerythritol tetranitrate single crystals: Mechanism of an isotropic sensitivity", Journal of Physical Chemistry A 104 11169 (2000). [102] V. K. Jindal, and D. D. Dlott, "Ori entation dependence of shockinduced heating in anharmonic molecular crystals ", Journal of Applied Physics 83 5203 (1998). [103] M. M. Kuklja, "On the initiation of ch emical reactions by electronic excitations in molecular solids", Applied Physics aMaterials Science & Processing 76 359 (2003). [104] M. M. Kuklja et al. "Role of electronic excitations in explos ive decomposition of solids", Journal of Applied Physics 89 4156 (2001). [105] M. M. Kuklja, and A. B. Kunz, "Com pressioninduced eff ect on the electronic structure of cyclotrimethylene trinitramine c ontaining an edge dislocation", Journal of Applied Physics 87 2215 (2000). PAGE 137 124 [106] M. M. Kuklja, and A. B. Kunz, "Ele ctronic structure of molecular crystals containing edge dislocations", Journal of Applied Physics 89 4962 (2001). [107] M. M. Kuklja, E. V. Stefanovich, and A. B. Kunz, "An excitonic mechanism of detonation initiation in explosives ", Journal of Chemical Physics 112 3417 (2000). [108] E. J. Reed et al. "A transient semimetallic laye r in detonating nitromethane", Nature Physics 4 72 (2008). [109] F. P. Bowden, and A. D. Yoffe, "Explosion in Liquids and Solids", Endeavour 21 125 (1962). [110] R. Menikoff, and T. D. Sewell, "Constituent properties of HMX needed for mesoscale simulations", Com bustion Theory and Modelling 6 103 (2002). [111] M. W. Conroy et al. "Firstprinciples anisotropic constitutive relationships in betacyclotetramethylene tetranitramine (b etaHMX)", Journal of Applied Physics 104 (2008). [112] M. W. Conroy et al. "Firstprinciples investigati on of anisotropic constitutive relationships in pentaerythritol te tranitrate", Physical Review B 77 (2008). [113] M. W. Conroy et al. "Density functional theory calculations of anisotropic constitutive relationships in alphacyclotrimethylenetrinitramine", Journal of Applied Physics 104 (2008). PAGE 138 ABOUT THE AUTHOR Michael W. Conroy earned a M.S. and B.S. in Physics at the University of South Florida (USF) following the receipt of his A.A. fr om Saint Petersburg Junior College. During his graduate studies, Michael wa s awarded a summer fellowship through the Naval Research Enterprise Internship Program by the American Society of Engineering Education at the Naval Research Laboratory (NRL). He also received three Fred L. and Helen M. Tharp Fellowships awarded by the De partment of Physics at USF. In addition, Michael has received several travel awards to present his research at scientific conferences from the American Physical Society Topical Gro up on Shock Compression of Condensed Matter and the USF Student Government. Following the completion of the Ph.D. degree, Michael began a Postdoctoral Fellowship at the NRL awarded by the National Research Council. xml version 1.0 encoding UTF8 standalone no record xmlns http:www.loc.govMARC21slim xmlns:xsi http:www.w3.org2001XMLSchemainstance xsi:schemaLocation http:www.loc.govstandardsmarcxmlschemaMARC21slim.xsd leader nam 22 Ka 4500 controlfield tag 007 crbnuuuuuu 008 s2009 flu s 000 0 eng d datafield ind1 8 ind2 024 subfield code a E14SFE0003196 035 (OCoLC) 040 FHM c FHM 049 FHMM 090 XX9999 (Online) 1 100 Conroy, Michael. 0 245 Density functional theory studies of energetic materials h [electronic resource] / by Michael Conroy. 260 [Tampa, Fla] : b University of South Florida, 2009. 500 Title from PDF of title page. Document formatted into pages; contains X pages. Includes vita. 502 Dissertation (Ph.D.)University of South Florida, 2009. 504 Includes bibliographical references. 516 Text (Electronic dissertation) in PDF format. 538 Mode of access: World Wide Web. System requirements: World Wide Web browser and PDF reader. 3 520 ABSTRACT: Firstprinciples calculations employing density functional theory (DFT) were performed on the energetic materials PETN, HMX, RDX, nitromethane, and a recently discovered material, nitrate ester 1 (NEST1). The aims of the study were to accurately predict the isothermal equation of state for each material, improve the description of these molecular crystals in DFT by introducing a correction for dispersion interactions, and perform uniaxial compressions to investigate physical properties that might contribute to anisotropic sensitivity. For each system, hydrostaticcompression simulations were performed. Important properties calculated from the simulations such as the equilibrium structure, isothermal equation of state, and bulk moduli were compared with available experimental data to assess the agreement of the calculation method. The largest contribution to the error was believed to be caused by a poor description of van der Waals (vdW) interactions within the DFT formalism. An empirical van der Waals correction to DFT was added to VASP to increase agreement with experiment. The average agreement of the calculated unitcell volumes for six energetic crystals improved from approximately 9% to 2%, and the isothermal EOS showed improvement for PETN, HMX, RDX, and nitromethane. A comparison was made between DFT results with and without the vdW correction to identify possible advantages and limitations. Uniaxial compressions perpendicular to seven lowindex crystallographic planes were performed on PETN, HMX, RDX, nitromethane, and NEST1. The principal stresses, shear stresses, and band gaps for each direction were compared with available experimental information on shockinduced sensitivity to determine possible correlations between physical properties and sensitivity. The results for PETN, the only system for which the anisotropic sensitivity has been thoroughly investigated by experiment, indicated a possible correlation between maximum shear stress and sensitivity. The uniaxial compressions that corresponded to the greatest maximum shear stresses in HMX, RDX, solid nitromethane, and NEST1 were identified and predicted as directions with possibly greater sensitivity. Experimental data is anticipated for comparison with the predictions. 590 Advisor: Ivan I. Oleynik, Ph.D. 653 Equation of state First principles Shear stress Explosives Compression 690 Dissertations, Academic z USF x Physics Doctoral. 773 t USF Electronic Theses and Dissertations. 4 856 u http://digital.lib.usf.edu/?e14.3196 