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

Full Text 
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 Ka controlfield tag 001 001935136 003 fts 005 20080421142313.0 006 med 007 cr mnuuuuuu 008 080421s2007 flua sbm 000 0 eng d datafield ind1 8 ind2 024 subfield code a E14SFE0002276 035 (OCoLC)225865327 040 FHM c FHM 1 100 Conroy, Michael W. 0 245 Firstprinciples studies of energetic materials h [electronic resource] / by Michael W. Conroy. 260 [Tampa, Fla.] : b University of South Florida, 2007. 520 ABSTRACT: Firstprinciples density functional theory calculations were performed on a number of important energetic molecular crystals, pentaerythritol tetranitrate (PETN), cyclotetramethylene tetranitramine (HMX), cyclotrimethylene trinitramine (RDX), and nitromethane. Simulations of hydrostatic and uniaxial compressions, as well as predictions of groundstate structures at ambient conditions, were performed using the DFT codes CASTEP and VASP. The first calculations done with CASTEP using GGAPW yielded reasonable agreement with experiment for the calculated isothermal EOS for PETNI from hydrostatic compression data, yet the EOS for betaHMX shows substantial deviation from experiment. Interesting anisotropic behavior of the shearstress maxima were exhibited by both crystals upon uniaxial compression.It was predicted that the [100] direction, the least sensitive direction of PETN, has significantly different values for shear stress maxima tauyx and tauzx, in contrast to the more sensitive directions, [110] and [001]. In addition, nonmonotonic dependence of one of the shear stresses as a function of strain was observed upon compression of PETN in the [100] direction. VASP calculations were later performed, and the results yielded good qualitative agreement with available experimental data for the calculated isothermal EOS and equilibrium structures for PETNI, betaHMX, alphaRDX, and nitromethane. Using VASP, uniaxial compression simulations were performed in the [100], [010], [001], [110], [101], [011], and [111] directions for all crystals up to the compression ratio V/V0 =0.70.The VASP calculations of PETN reproduced the CASTEP results of significantly different values of tauyx and tauzx for the insensitive [100] compression, and relatively high and equal values of tauyx and tauzx for the sensitive [110] and [001] compressions. A correlation between this behavior of shear stress upon uniaxial compression and sensitivity was suggested, and predictions of anisotropic sensitivity of HMX, RDX, and nitromethane were made. Further analysis of the VASP results for PETN do not indicate a correlation between sensitivity and shear stress maxima as a function of longitudinal stress, where longitudinal stress is an appropriate experimental independent variable for comparison. The validity of a correlation between shear stress maxima and sensitivity requires further investigation. Further characterization of the anisotropic constitutive relationships in PETN was performed. 502 Thesis (M.S.)University of South Florida, 2007. 504 Includes bibliographical references. 516 Text (Electronic thesis) in PDF format. 538 System requirements: World Wide Web browser and PDF reader. Mode of access: World Wide Web. 500 Title from PDF of title page. Document formatted into pages; contains 51 pages. 590 Advisor: Ivan I. Oleynik, Ph.D. 653 Materials modeling. Density functional theory. Equation of state. PETN. HMX. RDX. Nitromethane. 690 Dissertations, Academic z USF x Physics Masters. 049 FHMM 090 QC21.2 (ONLINE) 773 t USF Electronic Theses and Dissertations. 4 856 u http://digital.lib.usf.edu/?e14.2276 PAGE 1 FirstPrinciples Studies of Energetic Materials by Michael W. Conroy A thesis submitted in partial fulfillment of the requirements for the degree of Master of Science Department of Physics College of Arts and Sciences University of South Florida Major Professor: Ivan I. Oleynik, Ph.D. Dale Johnson, Ph.D. Lilia Woods, Ph.D. Date of Approval: October 26, 2007 Keywords: materials modeling, density function al theory, equation of state, petn, hmx, rdx, nitromethane Copyright 2007 Michael W. Conroy PAGE 2 Table of Contents List of Tables................................................................................................................. ....iii List of Figures......................................................................................................................v Abstract ...................................................................................................................... .......ix I. Introduction .............................................................................................................1 II. FirstPrinciples Calculations ...................................................................................3 A. A system of interacting electrons and nuclei ..............................................3 1. Kinetic energy of the nuclei ............................................................6 2. Electronelectron interactions .........................................................7 3. Electronnuclei interactions ............................................................8 B. Periodic systems ..........................................................................................9 C. Calculation of forces and stress.................................................................10 III. Summary of Results Obtained ..............................................................................11 i PAGE 3 IV. Energetic Materials at High Compression: FirstPrinciples Density Functional Theory and Reac tive Force Field Studies ...........................................13 V. Energetic Materials at High Compression: FirstPrinciples Density Functional Theory Studies ....................................................................................18 VI. Anisotropic Constitutive Relationships in Energetic Materials: PETN and HMX ...................................................................................................................27 VII. Anisotropic Constitutive Relationships in Energetic Materials: Nitromethane and RDX ........................................................................................32 VIII. FirstPrinciples Investigation of Anisotropic Constitutive Relationships in PETN ....................................................................................................................37 References .................................................................................................................... .....49 Appendix............................................................................................................................51 ii PAGE 4 List of Tables IV. Energetic Materials at Hi gh Compression: FirstPrinciples Density Functional Theory and Reactive Force Field Studies Table 1. Equilibrium lattice constants of PETN at zero temperature calculated by ReaxFF and DFT codes Castep and SeqQuest using several density functionals (LDA, GGA, GGAPBE) and compared with experiment..............................................................................................................16 V. Energetic Materials at High Compression: FirstPrinciples Density Functional Theory Studies Table 1. Equilibrium lattice parameters of PETNI calculated by DFT using several density functionals (LDA, GGAPW, and GGAPBE) and compared with experiment.....................................................................................21 iii PAGE 5 Table 2. Equilibrium lattice parameters of HMX calculated by DFT using several density functionals (LDA, GGAPW, and GGAPBE) and compared with experiment.....................................................................................21 VI. Anisotropic Constitutive Relationships in Energetic Materials: PETN and HMX Table 1. Elastic constants of PETNI cal culated in this work and compared with experiment.....................................................................................................30 VII. Anisotropic Constitutive Relationships in Energetic Materials: Nitromethane and RDX Table 1. Calculated and experi mental elastic constants of RDX..........................35 VIII. FirstPrinciples Investigation of An isotropic Constitutive Relationships in PETN Table 1. Calculated equilibrium structure of PETN compared with experiment..............................................................................................................41 Table 2. Bulk modulus results...............................................................................43 Table 3. Elastic constants of PETN.......................................................................43 iv PAGE 6 List of Figures IV. Energetic Materials at Hi gh Compression: FirstPrinciples Density Functional Theory and Reactive Force Field Studies Figure1. Crystal structure of PETNI (space group symmetry cP124 )................15 Figure 2. Constant pressure calc ulations of EOS for PETNI...............................17 Figure 3. Energy changes relative to the equilibrium zeropressure structure of PETNI upon uniaxial co mpression in [100] direction......................17 V. Energetic Materials at High Compression: FirstPrinciples Density Functional Theory Studies Figure 1. Isothermal hydrostatic equation of state of PETNI...............................22 Figure 2. Lattice parameters of PETNI as a function of compression ratio 0VV ......................................................................................................................22 v PAGE 7 Figure 3. Isothermal hydrosta tic equation of state of HMX...............................23 Figure 4. Lattice parameters of HMX as a function of compression ratio 0VV ......................................................................................................................23 Figure 5. Shear stresses in P ETNI upon uniaxial compression along [100], [110], and [ 001] directions as a func tion of compression ratio 0aa .........24 Figure 6. Shear stresses in HMX upon uniaxial compression along [110], [101], and [ 011] directions as a func tion of compression ratio 0aa .........25 VI. Anisotropic Constitutive Relationships in Energetic Materials: PETN and HMX Figure 1. Hydrostatic EOS of PETNI...................................................................29 Figure 2. Shear stresses xy and xz at uniaxial compression 0.7 in PETNI..........30 Figure 3. Hydrostatic EOS of HMX...................................................................30 Figure 4. Shear stresses xy and xz at uniaxial compression 0.7 in HMX..........31 vi PAGE 8 VII. Anisotropic Constitutive Relationships in Energetic Materials: Nitromethane and RDX Figure 1. Hydrostatic EOS of nitromethane..........................................................34 Figure 2. Shear stresses xy and xz at uniaxial compression 0.7 in nitromethane..........................................................................................................35 Figure 3. Hydrostatic EOS of RDX....................................................................35 Figure 4. Shear stresses xy and xz at uniaxial compression 0.7 in RDX...........36 VIII. FirstPrinciples Investigation of An isotropic Constitutive Relationships in PETN Figure 1. Isothermal hydrostatic EOS of PETNI..................................................42 Figure 2. Bond lengths of PETN under hydrostatic compression.........................42 Figure 3. Lattice parameters of P ETN under hydrostatic compression.................42 Figure 4. Stress tensor components xx yy and zz of PETN under uniaxial compression.............................................................................................44 vii PAGE 9 Figure 5. Change in energy per at om upon uniaxial and hydrostatic compression...........................................................................................................44 Figure 6. Band gap upon uniaxial and hydrostatic compression...........................44 Figure 7. Shear stress maxima xy and xz ..............................................................45 Figure 8. Maximum shear stress xy and xz as a function of longitudinal stress xx .................................................................................................................46 viii PAGE 10 FirstPrinciples Studies of Energetic Materials Michael W. Conroy ABSTRACT Firstprinciples density functional theory calculat ions were performed on a number of important energeti c molecular crystals, pentaerythritol tetranitrate (PETN), cyclotetramethylene tetranitramine (HMX), cyclotrimethylene trinitramine (RDX), and nitromethane. Simulations of hydrostatic and uniaxial compressions, as well as predictions of groundstate stru ctures at ambient conditions, were performed using the DFT codes CASTEP and VASP. The first calculations done with CA STEP using GGAPW yielded reasonable agreement with experiment for the calculated isothermal EOS for PETNI from hydrostatic compression da ta, yet the EOS for HMX shows substantial deviation from experiment. Interesting anisotropic behavior of the shearstress ma xima were exhibited by both crystals upon uniaxial compression. It wa s predicted that the <100> direction, the least sensitive direction of PETN, has signi ficantly different values for shear stress maxima yx and zx in contrast to the more sensitive directions, <110> and <001>. In addition, nonmonotonic dependence of one of th e shear stresses as a function of strain was observed upon compression of PETN in the <100> direction. ix PAGE 11 VASP calculations were later performed, and the results yielded good qualitative agreement with available experimental da ta for the calculated isothermal EOS and equilibrium structures for PETNI, HMX, RDX, and nitromethane. Using VASP, uniaxial compression simulations were perf ormed in the <100>, <010>, <001>, <110>, <101>, <011>, and <111> directions for all crystals up to the compression ratio V/V 0 = 0.70. The VASP calculations of PETN reproduced the CASTEP results of significantly different values of yx and zx for the insensitive <100> compression, and relatively high and equal values of yx and zx for the sensitive <110> and <001> compressions. A correlation between this behavior of sh ear stress upon uniaxial compression and sensitivity was suggested, and predictions of anisotropic sensitivity of HMX, RDX, and nitromethane were made. Further analysis of the VASP results for PETN do not indicate a correlation between sensitivity and shear stress maxima as a function of longitudinal stress, where longitudinal stress is an appropriate experi mental independent variable for comparison. The validity of a correlation between shear stress maxima and sensitivity requires further investigation. Further characte rization of the anisotropic c onstitutive relationships in PETN was performed. x PAGE 12 I. INTRODUCTION First principles calculations treat sy stems on the atomic level using quantum mechanics, with the intent of using as fe w approximations as feasibly possible. An exceptional theoretical tool for the implementation of firstprinciples calculations is densityfunctional theory (DFT), which th e development of DFT earned Walter Kohn a Nobel Prize in Chemistry in 1998. DFT was a re volutionary development in that it made accurate calculations possible for large syst ems, such as solids and large molecules. Numerous additional efforts to make calculations of this type possible have been undertaken, and the advance in computer hardware technology has ma de very accurate firstprinciples calculations on large syst ems possible in the last two decades. In particular, energetic materials, with large unit cells and complex molecular geometries, have only recently become possible to study via the firstprinciples methodology. The use of firstprinciples calculations ha s allowed for exceptional insight into the fundamental physics and chemistry involved in energetic materials research. Since experimentation with these materials can be difficult and costly, the use of theoretical simulation has proven to be a timeand cost effective strategy for the investigation of energetic materials. Importantly, the physical processes involved in the initiation of detonation in these materials at the atomic level are not we ll understood. Simulation from 1 PAGE 13 firstprinciples and comparison with expe rimental data may help elucidate the microscopic mechanisms responsible for the ex plosive reactions that these materials are known for. PETN (pentaerythritol tetranitrate) exhib its anisotropic shockinitiation sensitivity upon uniaxial compression, which was discovered by J. Dick 1 In fact, detonation has occurred when the crystal is compressed in sensitive directions at longitudinal stresses that may be simulated by firstprinciples ca lculations. By investig ating the calculated atomicscale properties upon compression in the se nsitive directions in contrast to the properties observed from compression in less sensitive directions, it may be possible to identify fundamental factors du e to intrinsic structural and chemical properties of PETN contributing to this strong anisotropy. Importa ntly, once mechanisms of sensitivity to detonation are identified, it is quite possible that these mechanisms may be extended to predict the initiation behavior in several othe r important energetic materials with similar properties, such as HMX, RDX, and nitromethane. This work presents firstprinciples calc ulations of the energetic materials PETN, HMX, RDX, and nitromethane. In Section II, a brief overview of the important ideas and methods involved in firstprinciples calcula tions is provided. Section III provides a summary of the motivation and goals of this work, as well as the important accomplishments and conclusions. Sections IVV III each contain articles that either have been published or will be submitted for publishing, and these articles summarize the work that we have done in the applic ation of firstprinciples simu lation to energetic materials. 2 PAGE 14 II. FIRSTPRINCIPLES CALCULATIONS The goal of employing firstprinciples me thods is to calculate the physical and chemical properties of materials by usi ng the fundamental equations of quantum mechanics with minimal approximations. Calc ulations of this type are also called abinitio calculations, meaning from the beginning in Latin. The purpose of this section is to briefly describe some of the important eq uations and approximations used in the firstprinciples calculations described in Sections IVVIII. A. A system of interacti ng electrons and nuclei The fundamental particles considered in a treatment of solids are the nuclei (neutrons and protons) and the electrons. A system of this type may be described by the manybody Hamiltonian 2 JI JI JI ji ji Ii Ii I I I I i i eeZZ e eZ M m H RR rr Rr2 2 2 2 2 2 22 1 2 1 2 2 (2.1) where the lowercase subscripts represent electrons and the uppercase subscripts are used for the nuclei. The first two terms repr esent the kinetic energy of the electrons and the nuclei, respectively. The third term represen ts the interaction of the electrons with the 3 PAGE 15 nuclei, the fourth is due to the electronelect ron interaction, and the fifth term describes the interactions between nuclei. A very brief description of the approximations used to handle this Hamiltonian is given below. 1. Kinetic energy of the nuclei The only term in (2.1) that can be consider ed small is the term that provides the kinetic energy of the nuclei 2 Since the mass of the nucleus is very large compared to the electrons, the quantity M 1 I is very small. Using the BornOppenheimer approximation, one may essentially freeze th e positions of the nuclei and abate the manybody problem by considering onl y the dynamics of the electrons 3 Hence, the kinetic energy of the nuclei would be neglected under th is approximation. Equation 2.1 can then be reduced to IIextEVVTH int (2.2) where denotes the kinetic en ergy of the electrons, describes the electronelectron interaction, is the interaction between the nuclei and the electrons (or the external potential), and is composed of the additional energy terms that are not inherent to the electrons, such as the classical nucleinuclei interaction T int V extV IIE 2 2. Electronelectron interactions 4 PAGE 16 The most difficult interaction to model is the electronelectron interaction. In the calculations that follow in Sections IVVIII, densityfunctional theory (DFT) was used to treat this interacti on. The two major works that formul ated DFT were by Hohenberg and Kohn 4 in 1964 and by Kohn and Sham 5 in 1965. In the first work, Hohenberg and Kohn 4 showed that there ex ists a unique external potential that determines a given groundstate charge density of a system of interacting electrons in the presence of nuclei, which was named the HohenbergKohn theorem. An important consequence of this theorem is that the groundstate charge density uniquely determines the groundstate energy of su ch a system, and hence the groundstate energy is a functional of the 3dime nsional chargedensity function r nEr n r V r n ext extV VT VVTnE int intr rrrr dVnnFext (2.3) where is a functional of r nF r n 6 describing the electronelectron interaction and the integral of the product of density and the ex ternal potential depict s the electronnuclei interaction. Approximations used in firstprinciples calculatio ns for dealing with the latter will be covered in the next section. When the functional r nE is minimized with respect to this yields the groundstate energy of a system. The advantage of this result is that the original manybody probl em is exactly reduced to finding the 3r n 5 PAGE 17 dimensional chargedensity th at minimizes the functional r nE However, r nF is not known, which does not allow the direct application of (2.3). In the second work, Kohn and Sham 5 showed that a system of interacting electrons with groundstate density r n can be mapped exactly to a system of noninteracting electrons with the same density, facilitating the solution of the problem. In this framework, singleelectron orb itals, named KohnSham orbitals j are the solutions of the selfconsistent KohnSham equations 6 rrrrrjj j xc H extsVVVT (2.4) and the energy functional (2.3) can be rewritten as 6 rrrrrrr dVnnEnEnTnEext xc H s (2.5) T s is the kinetic energy of the noninteracting electrons 6 j j j e sd m nT rrr r2* 22 2 (2.6) HE which includes electrostatic interacti ons, is known as the Hartree energy, rr rr rr r dd nne nEH22 (2.7) and which is known as the exchangecorrela tion energy, is made up of all the terms that remain xcE 6 In other words, the terms which provi de a large contribution to the total energy and are easy to calculate have been clearly expressed, and the rest of the terms are hidden in whose exact functional form is not known. Fortunately, there are approximations to the exchangecorrelation en ergy functional, such as the localdensity xcE 6 PAGE 18 approximation 7 (LDA) and the generalizedgradient approximation (GGA) with the PW91 8, 9 and PBE 10, 11 functionals, that have been su ccessful in describing physical properties of systems, such as energy diffe rences between structures and structural parameters, to within a few percent of experiment 3 In a firstprinciples calculation, the KohnSham equations are solved selfconsistently, such that an initial guess for th e density is made to calculate the potential used to solve the KohnSham equations, and th e resulting density is compared with the original guess for selfconsistency 2 If the two densities are not equal within a specified tolerance, the density is strategically update d and another cycle begins, and this whole process continues until selfconsistency is reached 2 3. Electronnuclei interactions The valence electrons participate in almost all of the important chemical and solidstate phenomena for which firstprinciples calculations are used to predict, but the core electron states stay almost unchanged. The use of pseudopotential theory 12 in firstprinciples calculations replaces the real pot ential experienced by the valence electrons from the nucleus and core states by a w eaker potential, called a pseudopotential. The weakening of the potential allows the wavefu nctions to be expanded in a much smaller set of planewaves, making the calculation of this interaction computationally feasible 3 Meanwhile, several of the important features of the real potential are maintained by the pseudopotential, and these vary depending on the type of pseudopotential used for a calculation. 7 PAGE 19 B. Periodic systems For periodic solids, the Bloch theorem st ates that the wavefunction for each electron can be written as the product of a periodic function and a plane wave 3 kexpjjfi rrk r (2.8) where the periodic function k jf r has the same periodicity as the underlying crystal lattice. The periodic por tion of the wavefunction k jf r may be expanded with a discrete set of plane waves as basis functions 3 k kexpjjfci GrG Gr (2.9) where G are the reciprocal lattice vectors of the crystal. Thus, the wavefunctions can be written as a sum of plane waves 3 G GkrGk r icj jexp, (2.10) This allows the electronic wavefunctions at each k point to be expanded in a discrete set of plane waves. Though discrete, the basis set would still need to be infinite to completely expand the wavefunc tion. However, the constants corresponding to lower kinetic energy are more important than those at higher kinetic energy Gk jc 3 A cutoff kinetic energy can be specified to truncate the basis set, allowing a finite basis to expand the electronic wave functions at each k point. In practice, one can determine an 8 PAGE 20 appropriate cutoff by increasing the kinetic en ergy cutoff of the basis set until the total energy of the system converges. The occupied states at each k point contribute to the electronic potential, and there are an infinite number of k points to consider in an infinite solid 3 However, the electronic wave functions do not change much in a region around a given k point. Thus, a single k point can be used to represent a small region of k space. Procedures for finding a special set of k points that will yield accurate total en ergies have been devised, and in the calculations that follow, the MonkhorstPack 13 method was used. Briefly, this method sets up a grid in the Brillouin zone wherein the intersections of the lines forming the grid are the k points used for a calculation. For practi cal purposes, the number of grid lines, and hence the density of k points, can be increased until convergence of the total energy is achieved. C. Calculation of forces and stress The HellmannFeynman theorem can be used to find the force conjugate to any parameter in the Hamiltonian that is used to describe a system 2 For instance, the force conjugate to the position of a nucleus R I is given by 2 I IE R F (2.11) 9 PAGE 21 By designating the Hamiltonian H e to include all the term s involving electrons and E II to describe the classical interac tion of the nuclei with one another, the derivative above can be expressed from firstorder perturbation theory as 2 I II I e e I I e IE H H H E RR R R R (2.12) Since, in the ground st ate, energy is a minimum with respect to all vari ations in the parameters of the wave function, the terms with derivatives of vanish. Thus, the force conjugate can be found from taking the derivativ es of terms with explicit dependence on the parameter R I 2 given by 2 I II I e IE H E R R R (2.13) In the calculations of Sections IVVIII, the stress tensor is calcul ated, and stress is a generalized force that may be calculated using the HellmannFeynman theorem as well. The stress tensor element ij is given by 2 ij ijE V 1 (2.14) where V is volume, ij is strain, and i j are Cartesian indices. Approximations are made in firstprinciples calculations to reduce the complexity of the interactions of electrons and nuclei vi a the application of ps eudopotentials. Further, these approximations have made firstprin ciples calculations of systems containing hundreds of atoms computationally feasible and the results yielded have shown remarkable agreement with experiment. 10 PAGE 22 III. SUMMARY OF RESULTS OBTAINED The firstprinciples calculations describe d in Sections IVVIII were performed on the energetic materials PETNI, HMX, RDX, and solid nitromethane. These materials are solids composed of organic mol ecules arranged in a lattice, and are thus called molecular crystals. Each of these crystals are high explosives 14 used in mainly military applications. In fact, PETN and RDX played an important role in the explosive compositions employed in World War II 15 These materials detonate at high pressures, which makes them excellent systems fo r the study of microscopic mechanisms responsible for initiation of detonation. Presently, these mechanisms are not well understood. Further, PETN exhibits a strong shockinduced det onation anisotropy 1, 1618 upon uniaxial compression. For exampl e, detonation has taken place in the <110> direction under shock strengths of about 4 GPa 17 whereas the first report of detonation in the least sensitive direction, <100>, has been re ported at a pressure of about 31 GPa 18 .Under the hypothesis that mechanochemical reactions may arise in regions of great shear stress in these materials, the investigation of a co rrelation between shear stresses and shockinitiation sensitivity was performed. The articles in Sections IVVIII report the results of this analysis. 11 PAGE 23 The major accomplishments of the work that follows are The equilibrium structures, hydro static compression, and uniaxial compressions of PETNI and HMX along selected crystallographic directions were calculated using the DFT code CASTEP 19 The CASTEP results indicate nonmonot onic dependence of shear stress maxima on strain in PETN for the inse nsitive direction <100>, but not for the sensitive direction, <110>. A correla tion is suggested and further work was started to investigate this possibility. The equilibrium structures, hydrostati c compressions and expansions, and uniaxial compressions and expansions in the <100>, <010>, <001>, <110>, <101>, <011>, and <111> directions of PETNI, HMX, RDX, and solid nitromethane were calculated using the DFT code VASP 20, 21 with the PBE functional 10, 11 and PAW potentials 22, 23 The VASP results indicate that the shearstress maxima at compression ratio V/V 0 = 0.7 are greater for uniaxial compression of PETN in the sensitive directions than for the in sensitive directions. Based on this observation, a suggestion of possi ble anisotropic shockinitiation sensitivity was made for HMX, RDX, and solid nitromethane. Upon further investigation of the sh earstress maxima as a function of longitudinal stress and comparison w ith experimental data, no clear correlation between shearstress maxima and shockinitiation sensitivity is indicated by the uniaxialco mpression data of PETN. 12 PAGE 24 IV. ENERGETIC MATERIALS AT HIGH CO MPRESSION: FIRSTPRINCIPLES DENSITY FUNCTIONAL THEORY AND REACTIVE FORCE FIELD STUDIES The following article 24 was published in the proceedings of the 14 th American Physical Society Topical Conference on Shock Compression of Condensed Matter and it has been reproduced here with permissions fr om the American Institute of Physics and the authors of the article. 13 PAGE 25 ENERGETIC MATERIALS AT HIGH COMPRESSION: FIRSTPRINCIPLES DENSITY FU NCTIONAL THEORY AND REACTIVE FORCE FIELD STUDIES I.I. Oleynik 1 M. Conroy 1 S.V. Zybin 2 L. Zhang 2 A.C. van Duin 2 W.A. Goddard III 2 and C. T. White 3 1 Department of Physics, University of South Florida, Tampa, FL 33620 2 Materials and Process Simulation Center, California Institute of Technology, Pasadena, CA 91125 3 Naval Research Laboratory, Washington, DC 20375 Abstract. We report the results of a comparative study of pentaerythritol tetranitrate (PETN) at high compression using classical reactive in teratomic potential ReaxFF and firstprinciples density functional theory (DFT). Lattice parameters of PETN I, the ground state structure at ambient conditions, is obtained by ReaxFF and two different density functional methods (plane wave and LCAO pseudopotential methods) and compared with experiment. Calculated energetics and isothermal equation of state (EOS) upon hydrostatic compression obtained by DFT and ReaxFF are both in good agreement with available experimental data. Our calculations of the hydrostatic EOS at zero temp erature are extended to high pressures up to 50 GPa. The anisotropic characteristics of P ETN upon uniaxial compression were also calculated by both ReaxFF and DFT. Keywords: Energetic materials, PETN, interatomic potentials, molecular dynamics, EOS PACS: 62.50.+p, 82.40.Fp, 81.30.Hd, 46.40.Cd INTRODUCTION Prediction of properties of energetic materials using atomicscale simulations techniques is one of the challenging areas of energetic materials (EM) research. Molecular dynamics (MD) simulation of EM using classical reactive interatomic potentials is a powerful modeling tech nique that is capable of addressing subnanometer and subpicosecond length and time scales of shock compression and detonation phenomena. However, the results of computer simulations can only be as reliable as the ability of the interatomic potentials to properly describe a variety of chemical effects including bondbreaking and bondmaking. Recently, the reactive force field ReaxFF has been developed based on fitting of an abinitio database of HCNO chemistry and is currently being actively used for MD simulations of EM [13]. One of the important issues is the transferability of ReaxFF, i.e. its ability to describe a rich chemistry and physics of energetic materials compressed at high pressures and temperatures. The reliable experimental data mostly exist only for static volumetric compressions that produce the isothermal equation of states (EOS) in a limited range of pressures and temperatures. Although this information is useful to validate reactive interatomic potentials, it is very limited in nature. Additional atomicscale information is urgently needed to understand the fundamental mechanisms 14 PAGE 26 of shock compression and detonation of EM materials. The modeling/simulation is able to provide such information that is very difficult or sometimes impossible to obtain from experiment. Density functional theory (DFT) has been very successful in recent years in simulating and predicting properties of a wide spectrum of materials from firstprinciples, i.e. without adhoc parameters that are usually present in empirical and semiempirical methods. In many cases the accuracy of DFT in si mulating properties of condensed phases is within a few percent compared to experiment. However, systems with weak van der Waals interactions, such as energetic molecular crystals, comprise a real challenge for DFT, because various density fu nctionals, including the local density approximation (LDA) and the generalized gradient approximation (GGA) face difficulties in describing systems having very small overlap of electronic densities from atoms constituting the system. Nevertheless, the DFT should work reasonably well for energetic materials at high pressures (because of substantial overlap of electronic densities from atoms) and can be a very efficient tool for generating a database of firstprinciples data that can be used for validation and fitting of interatomic potentials. The purpose of this work is to perform a comparative study of static high pressure properties of pentaerythritol tetranitrate PETN using both density functional theory (DFT) and ReaxFF. PETN is one of the classical energetic materials that has been extensively studied in recent years [4]. In addition to being an important secondary explosive, PETN exhibits an interesting property of strong anisotropy in the response to shock initiation of detonation [4]. Therefore, we decided to use PETN as a test bed for comparative studies using ReaxFF and DFT. We are particularly interested in equilibrium properties of the PETN stable crystalline phase at ambient conditions (PETNI), see FIGURE 1, and its equation of state (EOS) in a wide range of applied pressures, see FIGURE 2. COMPUTATIONAL DETAILS The reactive force field (R eaxFF) is a bondorder dependent interatomic potential that includes covalent interactions via bondorders. The bondorder is calculated from the instantaneous interatomic distances that are continuously updated in the course of simulation, thus allowing creation and breaking of chemical bonds. In addition to covalent interactions, ReaxFF includes Coulomb and van der Waals interactions traditionally present in classical force fields. However, these terms that are calculated for each atomic pair are effectively screened at short interatomic distances, thus reducing their contributio n in the region where covalent interactions dominate. The ReaxFF was fitted to an extensive database of molecular and crystal structures of CHNO containing species including reaction pathways crystal structures and binding energy curves for all possible bonds including bond angle and torsional angle dependences [13]. Recen tly, ReaxFF has been successfully applied to study thermal decomposition of nitromethane [2] and RDX [3]. FIGURE 1. Crystal structure of PETNI (space group symmetry 1P42C ). Density functional calculations were performed by using to different basis sets: plane waves (PW) and linear combination of atomic orbitals (LCAO). In both cases, high quality pseudopotentials were employed to remove th e core electrons from calculations, thus making firstprinciples calculations feasible. We carefully studied the completeness of the plane wave basis set by optimizing energy cutoff to get the convergence in energies, forces and stresses better than 0.01 eV/atom, 0.05 eV/ and 0.1 GPa respectively. We used Vanderbilt ultrasoft pseudopotentials (USPP) with 500 eV energy cutoff, additional calculations were performed with 700 eV to check the 15 PAGE 27 Table 1. Equilibrium lattice constants of PETN at zero temperature calculated by ReaxFF and DFT codes Castep and SeqQuest using several density functionals (LDA, GGAPW, and GGAPBE) and compared with experiment. PETNI a=b, (error, %) c, (error, %) c/a (error, %) Cell Volume, Experiment 9.38 6.71 0.71 589.5 LDA (Castep) 8.891 (5.2%) 6.453 (3.8%) 0.726 (+2.3%) 510.1 (13.5%) GGAPW (Castep) 9.600 (+2.3%) 6.796 (+1.3%) 0.708 (0.3%) 626.3 (+6.2%) GGAPBE (Castep) 9.820 (+4.7%) 6.950 (+3.6%) 0.708 (0.3%) 670.2 (+13.7%) GGAPBE (SeqQuest) 9.702 (+3.4%) 6.889 (+2.7%) 0.710 ( 0.0% ) 648.5 (+10.0%) ReaxFF 9.427 (+0.5%) 6.989 (+4.2%) 0.741 (+4.4%) 621.1 (+5.4%) convergence. The Brillouin kpoint sampling with kpoint density of 0.08 1 was chosen to get the energies, forces and stre sses converged to values cited above. We used the plane wave code Castep [5]. LCAO DFT calculations have been performed using LCAO code SeqQuest that uses normconserving pseudopotentials (NCPP) and highquality contracted Gaussian basis sets [6]. It is a computationally efficient code that enables very largescale calculations using rather modest computational resources. For present the calculations, the basis sets of doublezeta plus polarization (DZP) quality have been used with kpoint sampling of the Brillouin zone with similar density as used in planewave calculations. RESULTS AND DISCUSSION The critical issue in em ploying DFT methods to study molecular crystals is the choice of the appropriate density functional. We investigated the accuracy of several density functionals including standard local density approximation (LDA) and several generalized gradient approximation (GGA) functionals including PerdewWang (PW) and PerdewBurkeErnzerhof (PBE) by calculating the equilibrium lattice parameters of PETNI. Table 1 compares the results of DFT calculations using ReaxFF, CASTEP and PWUSPP and LCAO codes with experiment. As expected, the LDA gives strong overbinding of weak van der Waals forces which results in lattice parameters 5% smaller than experimental values. The best results were obtained using GGAPW in case of PWUSPP (errors less than 2.3%) and GGAPBE (errors less than 3.4%) for LCAONCPP. ReaxFF also predicts lattice constants in good comparison with experiment. The hydrostatic equation of states (EOS) for PETNI was obtained by performing constant pressure calculations at zero temperature. Castep data were generated using variable cell optimization under the cons traint of a diagonal stress tensor with fixed values of diagonal matrix elements equal to a desirable pressure. The ReaxFF calculations were performed using damped constant pressure molecular dynamics which is equivalent to conjugate gradient minimization of the total energy. This method also allowed the optimization of the unit cell of the crystal under the constant pressure constraint. The space symmetry of the crystal structure has been removed in order to relax symmetry imposed constraints. The atomic coordinates were also simultaneously optimized to get zero forces on atoms. Isothermal EOS, i.e. the dependence of pressure on volume, for PETNI obtained from the constant pressure simulations by Castep and ReaxFF are shown in FIGURE 2 and compared with experimental data by Olinger et al [7]. Both Castep and ReaxFF data compare well with each other as well as with experiment. The pressure range covered by our simulations extends beyond the experimental pressures up to 50 GPa. The energy changes in eV/atom relative to the equilibrium zero pressure structure of PETNI are also predicted by ReaxFF in good agreement with DFT results. We also performed an investigation of uniaxial compression of PETN. It is wellknown that PETN shows interesting sensitivity properties, i.e. strong 16 PAGE 28 FIGURE 2. Constant pressure ca lculations of EOS for PETNI. Top panel: energy changes relative to the equilibrium zero pressure structure of PETNI. Bottom panel: the isotherm al EOS of PETNI. FIGURE 3. Energy changes relative to the equilibrium zeropressure structure of PETNI upon uniaxial compression in [100] direction. anisotropy in the response in shock initiation of detonation. The shock loading conditions are characterized by fast uniaxial compression of the crystals along specific crystallographic directions. Therefore, we decided to investigate the predictions of both DFT and ReaxFF for uniaxial loading of the PETNI. The calculations were performed by creating a supercell where the PETN crystal is rotated to be oriented along a particular direction coinciding with xaxis of the unit cell, resulting in a nontetragonal unit cell. Then, the xdimension of the unit cell was strained by an appropriate scaling of the lattice parameter We calculated the characteristics of uniaxiallyloaded PETNI crystal in [100], [110], [010], and [001] directions. Due to space limitations we show only [100] results in Figure 3. We found that the DFTSeqQuest and ReaxFF are in very good agreement and it is not surprising, because ReaxFF was fitted using a database generated by SeqQuest. However, Castep calculations substantially deviate from both ReaxFF and SeqQuest. One of the sources of this disagreement is different GGA functionals: Castep calculations were done with PW but SeqQuest with PBE functionals. This results in slightly different lattice constants for the equilibrium crystal structures of PETN, see TABLE 1, which may affect the curvature of the a 0/ Eaa curve. We are currently investiga ting the sources of this disagreement and will report the full results in a subsequent publication. ACKNOWLEDGEMENTS IIO is supported by NSFNIRT (ECS0404137) and AROMURI (W901 1NF0510266). Funding at Caltech was provided by ONR and AROMURI. CTW is supported by ONR directly and through Naval Research Laboratory. REFERENCES 1. van Duin A.C.T., Dasgupt a S., Lorant F. and Goddard III W.A., J. Phys. Chem. A 105, 9396 (2001). 2. Strachan A. et al, Phys. Rev. Lett. 91, 098301 (2003). 3. Strachan A. et al, J. Chem. Phys. 122, 054505 (2005). 4. Dick J.J., J. Appl. Phys. 76, 2726 (1994); 81, 601 (1997). 5. Segall M.D. et al, J. Phys. Cond. Matt. 14, 2717 (2002). 6. Schultz P.A., SeqQuest el ectronic structure code: http://dft.sandia.gov/ Quest/SeqQ_Home.htm 7. Olinger B., Halleck P.M., and Cady H.H., J. Chem. Phys. 62 4480 (1975). 17 PAGE 29 V. ENERGETIC MATERIALS AT HIGH CO MPRESSION: FIRSTPRINCIPLES DENSITY FUNCTIONAL THEORY STUDIES The following article 25 was published in the Proceedings of the 13 th International Detonation Symposium and it has been reproduced here. This article is not copyrighted. 18 PAGE 30 ENERGETIC MATERIALS AT HIGH COMPRESSION: FIRSTPRINCIPLES DENSITY FUNCTIONAL THEORY STUDIES Ivan I. Oleynik Mike Conroy Sergey V. Zybin f and Carter T. White g Department of Physics University of South Florida, Tampa, FL 33620 f Materials and Process Simulation Center, California Institute of Technology, Pasadena, CA 91125 g Naval Research Laboratory Washington, DC 203755320 Abstract. In this paper we present results of systematic investigation of two important energetic molecular crystals, pentaerythritol tetranitrate (PETN, C[CH 2 ONO 2 ] 4 ) and cyclotetramethylene te tranitramine (HMX, [CH 2 NNO 2 ] 4 ) at high compression using firstprinciples density functional theory. The lattice parameters of PETNI, and HMX, the ground state structures at ambient conditions, were obtained and compared with experiment. The isothermal equations of state have been obtained at a wide range of compressions. In addition to accurate simulation of isotropic hydrostatic compression we have performed a series of uniaxial compressions in order to simulate the anisotropic response of EM under conditions close to those achieved in shock compressed EMs. The isothermal hydrostatic equation of state for PETNI was predicted by DFT in a reasonable agreement with experiment. However, the hydrostatic EOS for HMX as calculated by DFT using GGAPW shows substantial deviation from experiment. The uniaxial compression of both PETNI and HMX show interesting anisotropic behavior of the shear stresses. We found that the [100] direction, the least sensitive direction in PETNI, has nonequal x y and x z shear stresses in contrast to the behavior of shear stresses seen from two other directions, [110] and [001], which exhibit greater sensitiv ity. In addition, for the [100] direction, one of the shear stresses, x z exhibits nonmonotonic dependence on the uniaxial strain. This may be the reason that the [100] direction is less sensitive to shock initiation than all other crystallographic directions. For HMX, the shear stresses exhibit nonmonotonic behavior for all three directions. INTRODUCTION An important problem in energetic materials (EMs) research is the understanding of the mechanisms of sensitivity of explosives upon impact and other stimuli 13 The ultimate goal is to aid in developing safe explosives that will minimize both transportation costs an d accidental damage to both persons and property. One of the important steps in developing predictive capabilities is the 19 PAGE 31 generation of a knowledge base of fundamental materials properties from first principles based upon underlying atomic scale structure. Special attention is being focused on obtaining accurate equations of state (EOS) for several important classes of EMs. Due to the dominant role of shockinduced phenomena in developing the initial stages of detonation 2 such properties should be extended beyond traditional isotropic dependence of pressure upon volume to include stressdependent relationships that describe the anisotropic materials response upon dynamical loading. In this paper we present results of firstprinciples density functional theory investigation of two important energetic molecular crystals, pentaerythritol tetranitrate (PETN, C[CH 2 ONO 2 ] 4 ) 4 and cyclotetramethylen e tetranitramine (HMX, [CH 2 NNO 2 ] 4 ) 5 These systems have been already studied both experimentally 69 and theoretically 1014 In particular, isothermal compression of PETN has been investigated by Olinger at al 4 in 1975, but accurate experimental results on isothermal hydrostatic compression of HMX have appeared relatively recently 5,6 Interestingly, Dick experimentally discovered the anisotropy in sensitivity to shock initiation present in PETN 7 Less sensitivity information is available for HMX energetic crystals 3 Hydrostatic compression of PETN and HMX has been studied theoretically using classical interatomic potentials 10 firstprinciples DFT 11 and HartreeFock methods 12,13 In particular, Sorescu, Rice and Thompson investigated both PETN and HMX using a specially developed classical intermolecular potential for nitramines 10 A recent study of hydrostatic compression of PETN has been performed using DFT with a Gaussian basis set within an allelectron implementation 11 PETN 12 and HMX 13 crystals have also been studied by the periodic HartreeFock method. Byrd, Scuseria and Chabalowski have recently performed DFT studies of equilibrium crystal structures of several energetic molecular crystals using the planewave pseudopotential code VASP 14 Their main goal was to evaluate the accuracy of DFT in predicting properties of EM molecular crystals. This work is concerned with systematic investigation of two classical EM molecular crystals, PETN and HMX, at high compression using firstprinciples density functional theory as implemented in the planewave pseudopotential framework 15 In contrast to previous theoretical studies, isothermal equations of state have been obtained at a wide range of compressions. More importantly, in addition to accurate simulation of isotropic hydrostatic compression we have performed a series of uniaxial compressions in order to simulate the anisotropic response of EM under conditions close to those achieved in shock compressed EMs. We discuss the importance of shear stresses developed upon uniaxial compression for understanding the anisotropic sensitivity of EMs. COMPUTATIONAL DETAILS Firstprinciples density functional calculations were performed by using CASTEP, a planewave pseudopotential code 15 High quality ultrasoft pseudopotentials were employed to remove the core electrons from calculations, thus making firstprinciples calculations feasible. The advantage of a planewave basis set over the linear combination of atomic orbitals (LCAO) basis is in absence of basis set superposition errors adherent to LCAO methods as well as its ability to control the errors due to the finite size of the basis set. We carefully studied the completeness of the plane wave basis set by optimizing the energy cutoff (a parameter controlling the size of the basis set) to get the convergence in energies, forces and stresses better than 0.01 eV/atom, 0.05 eV/ and 0.2 GPa respectively. This convergence was checked for both uncompressed, equilibrium structures as well as highly compressed crystals. By performing test calculations at several energy cutoffs (300, 400, 500, 700 and 1000 eV) we found that the optimal cutoff is 500 eV. The Brillouin kpoint sampling with a kpoint density of 0.08 1 at the highest compression was chosen to get the energies, forces and stresses converged to the values cited above. The MonkhorstPack kpoint grid was kept fixed during the compression simulations to achieve a smooth change of physical properties upon compression. It is well known that most EM crystals have numerous polymorphs that are exemplified by phase transitions occurring upon compression. In this study we focus our attention on PETNI and HMX, the lowest energy crystal structures at ambient conditions. The PETNI phase has a tetragonal 142 P c unit cell, containing two PETN 20 PAGE 32 molecules, with 58 atoms in total 4 The phase of HMX has a monoclinic structure with two molecules per unit cell, resulting in 56 atoms in total 12/ Pc 5 EQUILIBRIUM PROPERTIES OF PETN AND HMX MOLECULAR CRYSTALS The critical issue in employing DFT methods to study molecular crystals is the choice of the density functional. Density functional theory (DFT) has been very successful in recent years in simulating and predicting properties of a wide spectrum of materials fr om firstprinciples. However, systems with weak van der Waals interactions, such as energetic molecular crystals, comprise a real challenge for DFT because various density functionals, including the local density approximation (LDA) and the generalized gradient approximation (GGA), face difficulties in describing systems having very small overlap of electronic densities from atoms constituting the system. In order to assess the performance of different density functionals in describing properties of EM molecular crystals we have investigated the equilibrium properties of PETNI and HMX by employing different flavors of the generalized gradient approximation (GGA) such as the PerdewWang (PW) and PerdewBurkeErnzerhof (PBE) functionals and the traditional local density approximation (LDA). Our strategy was to choose the density functional that predicts the equilibrium crystal properties closest to experimental values. The lattice parameters (lattice constants and cell angles) were optimized simultaneously with the atomic coordinates to have both zero forces on atoms and zero components of the stress tensor within the tolerances of 0.05 eV/ and 0.2 GPa respectively. Tables 1 and 2 compare the DFT results obtained using LDA, GGAPW and GGAPBE with experimental values for PETNI and HMX, respectively. As expected the LDA gives strong overbinding of weak van der Waals forces which results in lattice parameters 5% smaller than experimental values for both PETN and HMX. For the PETNI crystal, the best results were obtained using GGAPW (errors less than 2.3%). The LDA calculation for HMX shows the closest match to experiment, yet the absolute value of the percent errors for volume and lattice parameters obtained by LDA, GGAPW, and GGAPBE only differ slightly. Because our goal was to analyze the two TABLE 1. Equilibrium lattice parameterss of PETN I calculated by DFT using several density functionals (LDA, GGAPW, and GGAP BE) and compared with experiment. PETNI a=b, (error, %) c, (error, %) c/a (error, %) Cell Volume, Experiment 9.38 6.71 0.71 589.5 LDA 8.891 (5.2%) 6.453 (3.8%) 0.726 (+2.3%) 510.1 (13.5%) GGAPW 9.600 (+2.3%) 6.796 (+1.3%) 0.708 (0.3%) 626.3 (+6.2%) GGAPBE 9.820 (+4.7%) 6.950 (+3.6%) 0.708 (0.3%) 670.2 (+13.7%) TABLE 2. Equilibrium lattice parameters of HMX calculated by DFT using several density functionals (LDA, GGAPW, and GGAP BE) and compared with experiment. HMX a, (error, %) b, (error, %) c, (error, %) Cell Volume, Experiment 6.54 11.05 8.70 90 124.3 519.4 LDA 6.375 (2.5%) 10.474 (5.2%) 8.375 (3.7%) 90 123.7 465.4 (10.4%) GGAPW 6.755 (+3.3%) 11.497 (+5.2%) 9.096 (+4.6%) 90 124.7 580.9 (+11.8%) GGAPBE 6.760 (+3.4%) 11.505 (+4.1%) 9.095 (+4.6%) 90 124.8 581.1 (+11.9%) 21 PAGE 33 FIGURE 1. Isothermal hydrostatic equation of state of PETNI. Left panel simulation domain, right panel experimental domain of compression ratios V. 0/ V FIGURE 2. Lattice parameters of PETNI as a function of compression ratio V. Left panel simulation domain, right panel experime ntal domain of co mpression ratios V. 0/ V V 10/1 0/ structures within the same computational framework, we decided to use GGAPW for further calculations of hydrostatic and uniaxial compressions. Obviously, the errors are due to wellknown problems of DFT to describe properly the van der Waals interactions. VVstructures within the same computational framework, we decided to use GGAPW for further calculations of hydrostatic and uniaxial compressions. Obviously, the errors are due to wellknown problems of DFT to describe properly the van der Waals interactions. VV HYDROSTATIC COMPRESSION HYDROSTATIC COMPRESSION The hydrostatic compression of PETN and HMX was simulated using a twostep relaxation process. Starting from the compression ratio The hydrostatic compression of PETN and HMX was simulated using a twostep relaxation process. Starting from the compression ratio 0/ the lattice constants were scaled by the appropriate compression ra tio, but the fractional atomic coordinates were taken from the structure corresponding to the relaxed cell obtained from the previous stage in the compression. During step one, the atomic coordinates were relaxed, keeping the unit cell fixed. The resulting pressure after relaxation was read and used to form a diagonal, hydrostatic stress tensor. At stage two, all degrees of freedom including atomic coordinates, lattice constants and cell angles were relaxed under the constraint of a fixed stress tensor formed after step 22 PAGE 34 FIGURE 3. Isothermal hydrostatic equation of state of HMX. Left panel simulation domain, right panel experimental domain of compression ratios V. 0/ V FIGURE 4. Lattice parameters of HMX as a function of compression ratio V. Left panel simulation domain, right panel experime ntal domain of co mpression ratios V. 0/ V V V0/ V00.4/1.0 VV 0/ one within a tolerance of 0.2 GPa. Such a procedure allows fast convergence of relaxation calculations due to gradual changes of several physical parameters of the system. The space symmetry of the crystal structure has also been removed in order to relax symmetry imposed constraints. The compression ratio V was varied from 1.00 to 0.40 with the step size of 0.025. one within a tolerance of 0.2 GPa. Such a procedure allows fast convergence of relaxation calculations due to gradual changes of several physical parameters of the system. The space symmetry of the crystal structure has also been removed in order to relax symmetry imposed constraints. The compression ratio V was varied from 1.00 to 0.40 with the step size of 0.025. 0/ We first discuss results for PETN. Isothermal EOS, i.e. the dependence of pressure on volume obtained by DFT (GGAPW) within the compression ratio interval We first discuss results for PETN. Isothermal EOS, i.e. the dependence of pressure on volume obtained by DFT (GGAPW) within the compression ratio interval 00.4/1.0 VV is shown in Figure 1 and compared with experimental data by Olinger et al 6 The pressure range covered by our simulations extends beyond the experimental pressures up to 150 GPa. Our results show very good agreement with experiment when plotted together within the simulation compression interval 00.4/1.0 VV However, a closer examination within the range of volume ratios present in the experiment shows less agreement. From the equilibrium structure to a volume ratio of about 0.9, the pressure of the simulated crystal is greater than experiment. Between the volume ratios 0.90 and 0.84, very good agreement is shown. For lower volumes, there is a greater discrepancy between the calculated and experimental pressures. It is worth noting that the curvature of the experimental isotherm at is less than the calculated isotherm, which indicates that the theoretical bulk modulus, 0V/V = 1.0 / BVdPdV seems to be greater than the experimental value. The lattice parameters and show very good agreement with experiment, see Figure 2. The a c 23 PAGE 35 0.70.80.91Uniaxial compression, a/a0 0 2 4 6Shear stress, GPa xy xz [100]shear stress 0.70.80.91Uniaxial compression, a/a0 0 2 4 6Shear stress, GPa xy[100]shear stress xz 0.70.80.91Uniaxial compression, V/V 0 0 2 4 6Shear stress, G Pa xy[110]shear stress xz 0.70.80.91Uniaxial compression, V/V 0 0 2 4 6Shear stress, G Pa xy[110]shear stress[001]shear stress xz 0.70.80.91Uniaxial compression, a/a0 0 2 4 6 S hea r st r ess, G P a yz[001]shear stress xz 0.70.80.91Uniaxial compression, a/a0 0 2 4 6 S hea r st r ess, G P a yz xz greatest error is observed for the c parameter, which is roughly 2%, but the error for is much smaller. Meanwhile, at compression ratios below what has been observed experimentally, the lattice constants show abrupt changes at V/V a 0 = 0.65 and also at V/V 0 = 0.35 which may indicate the possibility of polymorphic phase transitions of the PETN crystal. For HMX, the agreement between the calculated and experimental EOS is not so good. This can be easily seen in from the graphs showing both the simulation and experimental domains of volume ratios, see Figure 3. One of the sources of this disagreement is errors in the equilibrium lattice parameters obtained by DFT. The compression simulations begin from an equilibrium structure having higher volume than the experimental volume. Therefore, the initial volume V may affect the curvature of the ~10% 0 0/ EVV curve. As far as lattice constants as a function of compression ratio are concerned, the agreement between theory and experiment is good, see Figure 4. UNIAXIAL COMPRESSIONS OF PETN AND HMX CRYSTALS We also performed an investigation of uniaxial compression of PETN and HMX crystals. The uniaxially compressed state of the crystal is directly related to the state the crystal experiences upon shock loading. At sufficiently large shock wave intensities, the time scale associated with the initial process of shock compression is on the order of picoseconds. Therefore, the lattice rapidly transforms at the shock wave front to a uniaxially compressed state. By inve stigating the mechanical properties of a uniaxially compressed crystal including shear stresses we have the possibility to explore the underlying atomicscale mechanisms of anisotropic sensitivity of EMs. It is well known that PETN shows interesting sensitivity properties, i.e. a strong anisotropy in the response to shock initiation of detonation. Dick, in his classical experiments, has found that PETNI is less sensitive in the [100] direction but more sensitive in the [110] direction upon impact initiation 9 Less is experimentally known about the sensitivity of HMX. The shock loading conditions are characterized by fast uniaxial compression of the crystals along specific crystallographic directions. Therefore, we decided to investigate the uniaxial loading of PETN and HMX along different crystallographic directions with the aim to correlate the anisotropic characteristics with their sensitivity. The calculations were performed by creating a supercell where both PETN and HMX crystals were rotated to be oriented along a particular direction coinciding with the xaxis of the unit cell, resulting in a nontetragonal unit cell. Then, the xdimension of the unit cell was strained by an appropriate scaling of the lattice parameter a. We calculated the stress characteristics of uniaxiallyloaded PETNI in the [100], [110], and [001] directions, see Figure 5. The shear stresses FIGURE 5. Shear stresses in PETNI upon uniaxial compression along [100], [110], and [001] directions as a function of com p ression ratio 0 / aa 24 PAGE 36 0. 5 0.60. 7 0. 8 0. 9 1 0 2 4 6 8 10Longitudinal stress, G Pa [110] shear stressesSxy Sxz 0.50.60.70. 8 0.91 0 2 4 6 8 10 L ong i t u d i na l s t ress, G P a [101] shear stressesSxy Sxz 0.50. 6 0. 7 0.80. 9 1 0 2 4 6 8 1 0 Longitudinal st r ess, GPa [011] shear sressesSxy Sxz FIGURE 6. Shear stresses in HMX upon uniaxial compression along [110], [101], and [011] directions as a function of compression ratio 0/ aa are our particular interest because they are usually considered to be the driving forces of plastic deformations in crystals. Several interesting observations have appeared. In particular, we have found that for both [110] and [001] directions, the x y and x z shear stresses are almost equal for all the states of a uniaxially compressed crystal. In contrast, x y and x z are different for the least sensitive direction, [100]. Interestingly, x z shows nonmonotonic dependence as a function of compression ratio 0/ aa The nonmonotonic behavior of shear stresses has been observed in other classes of shockcompressed material. For example, we have recently found that covale ntly bonded materials such as diamond and silicon exhibit an anomalous elastic response upon shock wave propagation: the weaker wave has both elastic and plastic regions of material behind the shock wave front, but at larger shock intensities the plastic deformations disappear. We related such unusual behavior with the nonmonotonic behavior of the shear stresses: the crystal undergoes all the stages of compression quickly until it reaches the final compression corresponding to very lo w values of the shear stresses. Because the shear stresses are the driving forces of plastic deformations, the compressions corresponding to the lowest values of shear stresses may be characterized by th e absence of the plastic deformations. This might be the reason of the absence of the initiation for the [100] crystallographic direction at the same level of compression (shock wave intensity) as for crystals compressed in [110] and [001] directions. The shear stresses in the uniaxiallycompressed HMX crystal also exhibit anisotropic, nonmonotomic behavior for all three directions, [110], [101] and [011], see Figure 6. There is no experimental information on the sensitivity of single crystal HMX. Therefore, more work, both experimental and theoretical, is needed to understand the sensitivity properties of HMX. CONCLUSIONS We have studied the PETNI and HMX molecular crystals using firstprinciples DFT implemented in the planewave pseudopotential framework. The equilibrium properties of PETN and HMX were studied using different density functionals. It was found that for PETN, GGAPW gives the smallest error ~2 for lattice constants and for unit cell volume compared to experiment, whereas the lattice parameters of HMX exhibit larger errors for all the DFT functionals. % ~6% ~5% The isothermal hydrostatic equation of state for PETNI was predicted by DFT in a reasonable agreement with experiment. However, the hydrostatic EOS for HMX as calculated by DFT using GGAPW shows substantial deviation from experiment. An obvious reason for this discrepancy is the substantially larger equilibrium volume of the unit cell as calculated by DFT which results in different curvature of the 0/ EVV curve as compared to experiment at the same absolute volume of the crystal. 25 PAGE 37 The uniaxial compression of both PETNI and HMX show interesting anisotropic behavior of the shear stresses. In PETNI we found that the [100] direction, the least sensitive direction has nonequal x y and x z shear stresses in contrast to the behavior of shear stresses seen from two other directions, [110] and [001], which exhibit greater sensitivity. In addition, for the [100] direction, one of the shear stresses, x z exhibits nonmonotonic dependence on the uniaxial strain. This might be a primary reason why the [100] direction is less sensitive to shock initiation than all other crystallographic directions. For HMX, the shear stresses exhibit nonmonotonic behavior for all three directions. Because no experimental information on shock sensitivity is available for this crystal, more work, both theoretical and experimental, is needed to clarify the fundamental mechanisms of shock sensitivity in EMs. REFERENCES 1. For recent review see special issue of Materials Science and Technology, R. Armstrong and J. Knott (Eds.), 22 Issue 4 (2006). 2. R.W. Armstrong, and W.L. Elban, Mat. Sci. Tech. 22 381 (2006). 3. S.M. Walley, J.E. Field, and M.W. Greenaway, Mat. Sci. Tech. 22 402 (2006). 4. A.D. Booth, and F.J.. Llewellyn, J. Chem. Soc., 837 (1947). 5. C. S. Choi and H. P. Boutin, Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem. 26 1235 (1970). 6. B. Olinger, P.M. Halleck, and H.H. Cady, J. Chem. Phys. 62 4480 (1975). 7. C.S. Yoo, and H. Cynn, J. Chem. Phys. 111 10229 (1999). 8. J.C. Gump, and S.M. Peiris, J. Appl. Phys. 97 053513 (2005). 9. J. J. Dick, Appl. Phys. Lett. 44 859 (1984); 10. D.C. Sorescu, B.M. Rice, and D.L. Thompson, J. Phys. Chem. B 103 6783 (1999). 11. C.K. Gan, T.D. Sewell, and M. Challacombe, Phys. Rev. B 69 035116 (2004). 12. H.V. Brand, J. Phys. Chem. B 109 13668 (2005). 13. F.J. Zerilli, and M.M. Kuklja, J. Phys. Chem. A 110, 5173 (2006). 14. E.F.C. Byrd, G.E. Scusseria, and C.F. Chabalowski, J. Phys. Chem. B 108, 13100 (2004). 15. M.D. Segal et al, J. Phys. Cond. Matt. 14 2717 (2002). 26 PAGE 38 VI. ANISOTROPIC CONSTITUTIVE RELATIONSHIPS IN ENERGETIC MATERIALS: PETN AND HMX The following article 26 was accepted to be published in the proceedings of the 15 th American Physical Society Topical Conf erence on Shock Compression of Condensed Matter and it has been reproduced here with perm issions from the American Institute of Physics and the authors of the article. 27 PAGE 39 ANISOTROPIC CONSTITUTI VE RELATIONSHIPS IN ENERGETIC MATERIALS: PETN AND HMX M. Conroy 1 I.I. Oleynik 1 S.V. Zybin 2 and C. T. White 3 1 Department of Physics, University of South Florida, Tampa, FL 33620 2 Materials and Process Simulation Center, California Institute of Technology, Pasadena, CA 91125 3 Naval Research Laboratory, Washington, DC 20375 Abstract. This paper presents results of firstpr inciples density functional calculations of the equation of state (EOS) of PETNI and HMX. The isotropic EOS for hydrostatic compression has been extended to include uniaxial compressions in the [100], [010], [001], [110], [101], [011], and [111] directions up to compression ratio V/V 0 = 0.70. Equilibrium properties, including lattice parameters and elastic constants, as well as hydrostatic EOS are in good agreement with available experimental data. The shear stresses of uniaxially compressed PETNI and HMX have been evaluated and their behavior as a function of compression ratio has been used to make predictions of shock sensitivity of these EMs. A comparison of predicted sensitivities with available experimental data has also been performed. Keywords: Energetic materials, PETNI and HMX, densityfunctional theory, equation of state, shock sensitivity, shear stresses PACS: 31.15.Ar, 31.15.Ew, 31.70.Ks, 62.50.+p, 64.30.+t INTRODUCTION One of the goals of energetic materials (EMs) research is to obtain accurate equations of state (EOS) for several important classes of EMs. Although a number of theoretical and experimental studies have been performed to obtain hydrostatic EOS for several EMs 19 there is an urgent need to extend isotropic EOS to include a description of the materials response upon uniaxial compression. The uniaxially compressed state of the crystal is directly related to the state that a crysta l experiences upon shock loading. The shear stre sses upon uniaxial compression are of particular interest because they are usually considered to be the driving forces behind plastic deformations in crystals and the sources of mechanochemical reactions behind the shock wave front. Anisotropic response to plane shock wave loading was observed in the experiments on PETN by J. J. Dick 10,11 Detonation at relatively low input stresses and run distances was observed for shocks delivered in the [110] and [001] crystallographic directions, yet detonation was not observed for shocks delivered in the [100] and [101] directions. The goal of this work is to obtain highly accurate hydrostatic EOS for PETNI and HMX and extend them to include uniaxial 28 PAGE 40 compressions in a wide range of compression ratios. We discuss the behavior of the shear stresses upon uniaxial compression and their potential importance in explaining the anisotropic shock sensitivity of PETNI and HMX. COMPUTATIONAL DETAILS The density functional calculations were performed using the Vienna AbInitio Simulation Package (VASP) 12 Special attention was paid to ensure high accuracy in the calculations. Test calculations were done to determine the types of density functionals and pseudopotentials, the value of the energy cutoff, and the density of kpoint sampling within the Brillouin zone that are adequate to reproduce available experimental data on equilibrium crystal structures. From these tests, the Generalized Gradient Approximation (GGA) density functional theory (DFT) with the PerdewBurkeErnzerhof (PBE) functional, the soft projected augmented wavefunction (PAW) pseudopotential with a 400 eV nominal cutoff, and a fixed MonkhorstPack grid corresponding to an average of 0.08 1 kspacing at maximum compression have been chosen. An energy cutoff of 700 eV was used, which yielded convergence better than 0.4 GPa in pressure, 0.0015 eV/atom in energy, and 0.015 eV/ in forces The equilibrium structures of PETNI and HMX were calculated by relaxing all the degrees of freedom of the experimental structures 4,13 including lattice parameters and atomic coordinates using the quasiNewton method as implemented in VASP. For all relaxations in this study, the convergence fo r electronic steps was set to 10 6 eV, and relaxation continued until the maximum force on any atom within the unit cell was less than 0.03 eV/. FIGURE 1. Hydrostatic EOS of PETNI. Hydrostatic compression calculations were performed on single unit cells of PETNI and HMX in the range of V/V 0 = 1.06 to 0.60 in increments of 0.02, where V 0 is the volume at zero pressure. Uniaxial compressions of both PETN and HMX were performed in the [100], [010], [001], [110], [101], [011], and [111] directions. For each compression direction, the calculated equilibrium unit cell was rotated such that the xaxis corresponded to the direction of compression. Subsequently, the xcomponent of each lattice vector was scaled by 2% up to V/V 0 = 1.06 and down to V/V 0 = 0.70. RESULTS AND DISCUSSION PETNI. The calculated unit cell parameters were compared to experimental data from Olinger, Halleck, and Cady 4 According to their work, the lattice constants a = b = 9.386 and c = 6.715 The lattice constants from this study (and the percent difference from experiment) were 9.617 (+2.5%), 9.612 (+2.4%), and 6.83 (+1.7%) for a, b, and c, respectively. The disagreement between DFT calculations and experiment for molecular crystals has been studied in several works 2,3,14 and has been attributed to poor description of van der Waals interactions by DFT. The isothermal equation of state (EOS) from the hydrostatic compression simulations of this work was compared with the experimental data of Olinger et al 4 see Figure 1. Pressure, change in energy per atom, and band gap were plotted as a function of absolute volume in the range of V/V 0 = 1.0 to 0.60. Good qualitative agreement was observed between the EOS calculated in this study and experimental data. For PETNI, the elastic constants C 11 C 22 and C 33 were calculated and compared with experiment 15 see Table 1. Results on isotropic EOS and elastic constants 29 PAGE 41 TABLE 1. Elastic constants of PETNI calculated in this work and compared with experiment 15 C 33 (GPa) Study C FIGURE 2. Shear stresses yx and zx at uniaxial compression 0.7 in PETNI. FIGURE 3. Hydrostatic EOS of HMX. 11 (GPa) C 22 (GPa) Exp.: 17.22 17.22 12.17 Ref. 6 This work 18.25 18.5 14.5 (+6.0%) (+7.4%) (+19%) have given us confidence that DFT gives a reasonable description of the compressed crystals despite the problem with van der Waals interactions. In the case of uniaxial compressions of PETNI, the change in energy per atom, longitudinal stress, band gap, and shear stresses were examined as a function of c/c 0 The shear stresses, which showed the greatest relative variation among the inspected properties, were calculated from the equation 2xx yzx yyzz where x indicates the direction of compression. Figure 2 shows the values of the shear stresses at uniaxial compression of c/c 0 = 0.7. From the experiments of J. Dick 10 on PETNI, the [110] and [001] directions have exhibited a greater sensitivity to detonation via shock, and the [100] and [101] directions have displayed decreased sensitivity. It can be observed from Figure 2 that the shear stresses for the more sensitive directions, [110] and [001], exhibit greater values for both yx and zx than the directions of lower sensitivity, [100] and [101]. Reduced values of shear stresses and greater separation between yx and zx were observed in the less sensitive directions, [100] and [101]. Therefore, the suggestion 10 that shear stress behavior is linked with sensitivity is supported by these results. Nonmonotonic dependence of shear stress upon strain was also observed for directions of decreased sensitivity. From these calculations and experimental data, we see that greater calculated values of both yx and zx upon uniaxial compression correlates positively with increased sensitivity in energetic molecular crystals. HMX. The calculated equilibrium lattice constants of HMX were compared with the experimental data from Choi and Boutin 13 Their results were a = 6.54 b = 11.05 and c = 8.70 The calculated values of the lattice constants from this work (and the percent difference from experiment) were 6.701 (+2.5%), 11.347 (+2.7%), and 8.910 (+2.4%) for a, b, and c, respectively. Again, overestimation of the lattice parameters was associated with the inability of DFT functionals to describe the significant intermolecular van der Waals interactions. The calculated isothermal EOS was compared with the experimental data obtained by Gump and Peiris 5 and Yoo and Cynn 6 see Fig. 3. As observed in PETNI, DFT provides good qualitative agreement with the EOS data gathered 30 PAGE 42 FIGURE 4. Shear stresses yx and zx at uniaxial compression 0.7 in HMX. from experiment for HMX. The shear stresses for the uniaxial compressions performed on HMX at c/c 0 = 0.7 were plotted in Fig. 4. Greater magnitudes of shear stresses, which corre sponded to directions of greater sensitivity in PETNI, were present in the directions [110] and [010]. Consequently, if the trend found for PETNI is maintained for HMX, then greater sensitivity to initiation is expected in these directions. CONCLUSIONS The energetic material crystals PETNI and HMX were studied using firstprinciples density functional theory. Good agreement was observed between available experimental data and calculated equilibrium structures, hydrostatic equations of state, and elastic constants. The errors were attributed to poor description of van der Waals interactions by DFT. The inclusion of van der Waals interactions will be the subject of future work. The isotropic EOS have also been extended to include uniaxial compressions PETNI and HMX in the [100], [010], [001], [110], [101], [011], and [111]. Upon comparison of the simulated uniaxial compression data with experimental information on the anisotropy in sensitivity to shockinduced detonation for PETNI, a correlation between shear stresses and sensitivity was suggested. The magnitudes of the calculated shear stresses were greater for the more sensitive directions, [110] and [001], than for the less sensitive directions, [100] and [101]. Extrapolating this trend between calculated shear stresses and known sensitivity in PETN to our calculated results for HMX suggests that for HMX the [110] and [010] directions could be more sensitive to initiation, ACKNOWLEDGEMENTS The work at USF was supp orted by the Office of Naval Research (ONR) through the Naval Research Laboratory (NRL). The work at NRL was also supported by ONR both directly and through NRL. M. Conroy thanks the organizing committee of SCCM 2007 (Ricky Chau) for the travel support to attend SCCM 2007 in Hawaii. REFERENCES 1. Sorescu D. C., Rice B. M., and Thompson D. L., J. Phys. Chem. B 103, 6783 (1999). 2. Byrd E. F. C., and Rice B. M., J. Phys. Chem. C 111, 2787 (2007). 3. Liu H., Zhao J., Wei D., and Gong Z., J. Chem. Phys. 124 124501 (2006). 4. Olinger B., Halleck P.M. and Cady H.H., J. Chem. Phys. 62, 4480 (1975). 5. Gump J. C., and Peiris S. M., J. App. Phys. 97, 053513 (2005). 6. Yoo C., and Cynn H., J. Chem. Phys. 111, 10229 (1999). 7. Gan C., Sewell T. D., and Challacombe M., Phys. Rev. B 69, 035116 (2004). 8. Brand H. V., J. Phys. Chem. B 109 13668 (2005). 9. Sewell T. D., Menikoff R., Bedrov D., and Smith G. D., J. Chem. Phys. 119 7417 (2003). 10. Dick J. J., Appl. Phys. Lett. 44, 859 (1984). 11. Dick J.J., J. Appl. Phys. 76, 2726 (1994); 81, 601 (1997). 12. Kresse, G.; Furthmuller, J. Vienna AbInitio Simulation Package (VASP) 13. Choi C. S., and Boutin H. P., Acta Crystallogr., Sect. B 26 1235 (1970). 14. Byrd E. F. C., Scuseria G. E., and Chabalowski C. F., J. Phys. Chem. B 108, 13100 (2004). 15. Winey J. M., and Gupta Y. M., J. App. Phys. 90 1669 (2001). 31 PAGE 43 VII. ANISOTROPIC CONSTITUTIVE RELATIONSHIPS IN ENERGETIC MATERIALS: NITROMETHANE AND RDX The following article 27 was accepted to be published in the proceedings of the 15 th American Physical Society Topical Conf erence on Shock Compression of Condensed Matter and it has been reproduced here with perm issions from the American Institute of Physics and the authors of the article. 32 PAGE 44 ANISOTROPIC CONSTITUTI VE RELATIONSHIPS IN ENERGETIC MATERIALS: NI TROMETHANE AND RDX I.I. Oleynik 1 M. Conroy 1 and C. T. White 2 1 Department of Physics, University of South Florida, Tampa, FL 33620 2 Naval Research Laboratory, Washington, DC 20375 Abstract. The anisotropic constitutive relationships in solid nitromethane (NM) and RDX were studied using firstprinciples dens ity functional theory (DFT). In addition to hydrostatic compressions, we performed uniaxial compressions in the [100], [010], [001], [110], [101], [011], and [111] directions up to the compression ratio V/V 0 = 0.70. Equilibrium properties, including lattice parameters and elastic constants, as well as hydrostatic EOS, are in good agreement with available experimental data. The shear stresses of uniaxially compressed NM and RDX were used to predict the relative shock sensitivity between different crystallographic directions. Keywords: Energetic materials, nitromethane and RDX, densityfunctional theory, equation of state, shock sensitivity, shear stresses PACS: 31.15.Ar, 31.15.Ew, 31.70.Ks, 62.50.+p, 64.30.+t INTRODUCTION Shock sensitivity of energetic materials is one of the outstanding scientific issues that have been a focus of the energetic materials (EMs) community in the last se veral decades [14]. A number of experiments have investigated shock sensitivity of single crystals of PETNI [24] and HMX [5]. In particular, J.J. Dick observed the anisotropic response of PETN to plane shock wave loading: the [110] and [001] crystallographic directions were found to be more sensitive than the [100] and [101] directions [24]. Less is known about the relative shock sensitivity of other EMs such as solid nitromethane (NM) and RDX. Ultimately, the shock sensitivity must be related to the underlying atomicscale structure of EMs which governs their constitutive properties. The goal of this work is to provide such constitutive relationships for NM and RDX from firstprinciples by explicitly treating electrons using density functional theory (DFT). By extending the hydrostatic equations of state (EOS) to include a description of the materials response upon uniaxial compression, we are trying to establish the relationship between the behavior of the shear stresses and the anisotropic shock sensitivity of EMs. Such uniaxially compressed states of the crystal are directly related to the state that a crystal experiences upon shock loading. The shear stresses upon uniaxial compression are usually considered to be the driving forces behind plastic deformations in crystals and the source s of mechanochemical reactions behind the shock wave front. 33 PAGE 45 COMPUTATIONAL DETAILS FIGURE 1. Hydrostatic EOS of nitromethane. The firstprinciples calculations of EM crystals NM and RDX were performed using the density functional theory code VASP (Vienna AbInitio Simulation Package) [6]. Special attention was paid to ensure the accuracy of the calculations by optimizing parameters such as the types of density functionals and pseudopotentials, the energy cutoff and the density of kpoint sampling of the Brillouin zone. We found that the Generalized Gradient Approximation (GGA) DFT with the PerdewBurkeErnzerhof (PBE) functional, projected augmented wavefunction (PAW) pseudopotentials, 700 eV energy cutoff and fixed MonkhorstPack grid corresponding to an average of 0.08 1 kspacing at maximum compression are adequate to reproduce available experimental data on equilibrium properties of NM and RDX. The equilibrium stru ctures of NM and RDX were calculated by relaxing all degrees of freedom including lattice parameters and atomic coordinates using the quasiNewton method as implemented in VASP. For all relaxations in this study, the convergence fo r electronic steps was set to 10 6 eV, and relaxation continued until the maximum force on any atom within the unit cell was less than 0.03 eV/. Hydrostatic compression calculations were performed on single unit cells of NM and RDX in the range of V/V 0 = 1.06 to 0.70 in increments of 0.02. Uniaxial compressions were performed in the [100], [010], [001], [110], [101], [011], and [111] directions. For each compression direction, the calculated equilibrium unit cell was rotated such that the xaxis corresponded to the direction of compression. Subsequently, the xcomponent of each lattice vector was scaled by 2% up to V/V 0 = 1.06 and down to V/V 0 = 0.70. RESULTS AND DISCUSSION Nitromethane. The experimental equilibrium unit cell structure of solid nitromethane is orthorhombic, with space group P2 1 2 1 2 1 and lattice parameters , [7]. The calculated equilibrium lattice parameters are , which are within and 5.183 a 6.236 b 8.518 c 5.301 a 6.591 b 8.838 c 2.3% 5.7% 3.8% of experimental values, respectively. The source of these errors is in the poor description of van der Waals interactions by DFT. The isothermal equation of state (EOS) for the hydrostatic compression calculated by DFT was compared with the experimental work of Yarger, and Olinger [8], see Figure 1. Pressure, change in energy per atom, and band gap were plotted as a function of absolute volume in the range of V/V 0 = 1.0 to 0.70. Good qualitative agreement was observed between the EO S calculated in this study and experimental data. The anisotropic properties of NM have been investigated by performing uniaxial compressions in the [100], [010], [001], [110], [101], [011], and [111] directions up to compression ratio c/c 0 = 0.70. For each compression ratio, the change in energy per atom, longitudinal stress, band gap, and shear stresses were determined. The shear stresses were calculated as 2xx yzx yyzz where x x and () yyzz are the longitudinal and transverse stresses, respectively. The calculated values of yzx 34 PAGE 46 FIGURE 2. Shear stresses yx and zx at uniaxial compression 0.7 in nitromethane. FIGURE 3. Hydrostatic EOS of RDX. correspond to maximum sh ear stresses applied at to the compression direction. Figure 2 shows the values of the shear stresses at uniaxial compression c/c o45 0 = 0.7. The calculated isothermal EOS for RDX is shown in Figure 3. The only experimental data on the isothermal EOS of RDX were obtained by Olinger, Roof and Cady in 1978 [11] and only pressures up to 4 GPa were reported. Therefore, we do not plot these experimental points. By inspecting the values of the shear stresses along different crystallographic directions, we found that the [001] compression has the highest values of yx and z x In our work on uniaxial compressions of PETNI [9 ], we found a positive correlation between directions exhibiting greater sensitivity in shockco mpression experiments and the crystallographic directions possessing relatively higher shear stresses in theoretical compression simulations. Similar correlation between the directions of the smallest shear stresses and the less sensitive directions was found. Although there is no anisotropic shock sensitivity data available for NM, by extrapolating the correla tion found in PETNI, suggests that the [001 ] direction will be the direction of the highest shock sensitivity. An important method to test the accuracy of DFT is the comparison of predicted elastic constants with experiment. For RDX, the elastic constants C 11 C 22 and C 33 were calculated and compared with experiment [12], see Table 1. The average error of DFT is 3.2%, which gives us confidence that DFT pr ovides a sufficiently accurate description of the compressed state of EMs despite the problem with van der Waals interactions at small compression. The shear stresses for the uniaxial compressions of RDX at c/c 0 = 0.7 are plotted in Figure 4. Greater magn itudes of shear stresses, RDX. The experimental crystal structure of RDX possesses the orthorhombic space group Pbca with lattice parameters , [10]. The DFT results are , and the percent difference from experiment are +2.3%, +1.6%, and +3.7% respectively. Again, over estimation of the lattice parameters was associated with the inability of DFT functionals to describe intermolecular van der Waals interactions. 13.182a 11.574 b 10.709 c 13.522a 11.758b 11.102c TABLE 1. Calculated and experimental [12] elastic constants of RDX. Constants C 11 (GPa) C 22 (GPa) C 33 (GPa) This work 25.5 20.5 18.0 (0.4%) (3.8%) (5.3%) Exp.: 25.6 21.3 19.0 Ref. 12 35 PAGE 47 FIGURE 4. Shear stresses yx and zx at uniaxial compression 0.7 in RDX. which could correspond to directions of greater sensitivity, were present in the [100] and [011] directions of RDX. Thus, increas ed sensitivity for these directions could be expected. The [101] direction displayed the smallest shear stresses suggesting that it could be the least sensitivity. Again, these suggestions are made following the establishment of a positive correlation between experimental sensitivity data and shear stress values in PETNI [9]. CONCLUSIONS We have studied the energetic materials crystals NM and RDX using firstprinciples density functional theory and found good agreement between available experimental data and calculated equilibrium structures, hydrostatic equations of state, and elastic constants. The errors are probably due to an inadequate description of van der Waals interactions by DFT. The isotropic EOSs have been extended to include uniaxial compressions in the [100], [010], [001], [110], [101], [011], and [111]. Using a correlation found in previous studies of PETNI, we also made suggestions of relative shock sensitivities of NM and RDX based on calculated shear stresses for different crystallographic directions. For NM, the magnitudes of the shear stresses were the greatest for the [001] direction suggesting that it corresponds to the most se nsitive direction. In the case of RDX, the [100] and [011] directions are expected to be more se nsitive, while the [101] direction is expected to be the least sensitive direction. These theore tical predictions, which assume that the positiv e correlation between calculated shear st ress and observed sensitivity in PETN can be extrapolated to NM and RDX, require experimental confirmation. ACKNOWLEDGEMENTS The work at USF was supported by the Office of Naval Research (ONR) through the Naval Research Laboratory (NRL). The work at NRL was also supported by ONR both directly and through NRL. M. Conroy thanks the organizing committee of SCCM 2007 (Ricky Chau) for the travel support to attend SCCM 2007 in Hawaii. REFERENCES 1. Walley S.M., Field J.E., and Greenaway M.W., Mat. Sci. and Tech. 22 402 (2006). 2. Dick J.J., Mulford R.N., Spencer W. J., Pettit D.R., Garcia E., and Shaw D.C., J. Appl. Phys. 70, 3572 (1991). 3. Dick J.J. and Ritchie J.P., J. Appl. Phys. 76, 2726 (1994). 4. Dick J.J., J. Appl. Phys. 81, 601 (1997). 5. Hooks D.E., Ramos K.J., in Proceedings of the 13 th International Detonation Symposium, Norfolk, (2006). 6. Kresse, G.; Furthmuller, J. Vienna AbInitio Simulation Package (VASP). 7. Trevino S.F., Prince E., and Hubbard C.R., J. Chem. Phys. 73, 2996 (1980). 8. Yarger F.L., and Olinger B, J. Chem. Phys. 85 1534 (1986). 9. Conroy M., Oleynik I.I, Zybin S., and C.T. White, to be published in this volume. 10. Choi C.S., and Prince E., Acta Cryst. B 28 2857 (1972). 11. Olinger, B. W.; Roof, B.; Cady, H. H. Symposium international Sur Le Comportement Des Milieus Denses Sous Hautes Pressions Dynamiques; Commissariat a lEnergie Atomique Cnetre dEtudes de Vajours: Paris, France, 3 (1978). 12. Schwarz R.B., Hooks D.E., Dick J.J., Archuleta J.I., and Martinez A.R., J. Appl. Phys. 98, 056106 (2005). 36 PAGE 48 VIII. FIRSTPRINCIPLES INVESTIGATION OF ANISOTROPIC CONSTITUTIVE RELATIONSHIPS IN PETN The following article will be submitted to Physical Review B 37 PAGE 49 FIRSTPRINCIPLES INVESTIGATION OF ANISOTROPIC CONSTITUTIVE RELATIONSHIPS IN PENTAERYTHRITOL TETRANITRATE (PETN) M. Conroy I. I. Oleynik S. V. Zybin and C. T. White Department of Physics, Universi ty of South Florida, Tampa, FL 33620 Materials and Process Simulati on Center, California Institut e of Technology, Pasadena, CA 91125 Naval Research Laboratory, Washington, DC 20375 Abstract. Firstprinciples density functiona l theory (DFT) calculations have been used to obtain the c onstitutive relationships of the energetic material crystal, pentaerythritol tetranitrate (PETNI). The isotropic equation of state (EOS) for hydrostatic compression has b een extended to include uniaxial compressions in the <100>, <010>, < 001>, <110>, <101>, <011>, and <111> crystallographic directions up to a compression ratio of DFT predicts equilibrium properties such as lattice parameters and elastic constants, as well as hydrostatic EOS in agreement with available experimental data. The hydrostatic EOS has been extended to include the response of EMs upon anisotropic uniaxial compressions. Our resu lts show a substan tial anisotropy of various properties of PETNI upon uniax ial compressions. To characterize the anisotropic properties of PETN, different physical pr operties of the uniaxially compressed crystal such as the energy per atom, band gap, and stress tensor have been evaluated as a function of compression ratio. The maximum shear stresses were calculated and examined for a correlation with the anisotropy in shockinitiation sensitivity. 0/0. VV 7 I. INTRODUCTION One of the goals of energetic materials (EMs) research is to investigate the fundamental physical and chemical properties of EM molecular cr ystals from firstprinciples based on the underlying atomic scale structure. Accurate equations of state (EOS) are of particular interest because they establish fundamental relationships between thermodynamic variables and provide necessary input for the de scription of materials at the mesoscopic and continuum levels. Ultimately, th ese EOS are governed by the inte ractions at the atomic scale which provides an excellent opportunity for atom icscale simulation techniques to establish these relationships by linking the microscopic and macroscopic properties of the material. Most of the theoretical and experimental investigations focus on isotropic EOSs, i.e. relationships between pressure and volume. Ho wever, due to the dominant role of shockinduced phenomena during initiation and propagati on of detonation, there is an urgent need to extend the isotropic EOS to include the response of EMs upon anisotropic uniaxial compressions. PETN, an important EM, has attracted consid erable interest from the EM community because of the discovery of its anis otropic shock sensitivity by J.J. Dick 13 Some prior work 38 PAGE 50 on PETN has indicated anomalous behavior up on shock compression. In particular, Halleck and Wackerle 4 performed shock compressions of up to about 4 GPa on singlecrystal samples with planes cut perpendicular to the <110> and <001> directi ons. For shocks at and above 4 GPa on the <110> crystals, shockinduced deco mposition was expected from a rapid increase in stress amplitude approximately 0.3 s after shock arrival. However, it was J.J. Dick who unambiguously demonstrated in his classic e xperiments the anisotropy in shock sensitivity. He found that when compressed in its most sensitive direction, <110>, detonation has been observed 3 in PETN crystals under shock stre sses as small as 4.2 GPa (from J. Dick 3 < hkl > denotes the direction pe rpendicular to the ( hkl ) plane). The <001> direction has also been reported as sensitive 13 However, it is very insensitive to shocks in the <100> direction 13, 5 and the <101> direction has also been found to be relatively insensitive 1, 3 Later, the anisotropic shock sensitivity of PETN was also investigated by Soulard 6, 7 There have also been a number of works suggesting mechanisms for the directional dependence of sensitivity in PETN. Dick et al. 13 have proposed a model based upon steric hindrance to shear, which suggests that sensitiv ity is related to the ease by which slip is allowed via the available slip systems under shock compression. Further, Jindal and Dlott 8 have suggested that anisotropic heating in response to shock compression might be responsible for directional dependen ce of sensitivity. Gruzdkov and Gupta 8 have added to the possible mechanisms by arguing that changes in th e local polarization of the lattice arise due to conformational changes in PETN mol ecules under compression, allowing decomposition via ionic reactions. The macroscopic properties of several en ergetic materials in cluding PETN, HMX, RDX, and TATB have been studied and described well 9, 10 However, the microscopic behavior responsible for the shock initiation observed on the macroscopic scale is not well understood. Firstprinciples simula tion of these materials at the atomic and molecular level is therefore providing an excellent opportunity for the investigation of microscopic processes leading to shock initiation. Becau se of the relatively large un it cells and complex geometries of these energetic materials, firstprinciples ca lculations have only been feasible for the last decade due to both advances in computati onal methodology and the exponential growth in computational power of modern computers. PETN is an exempl ary EM for these calculations, as its anomalous anisotropi c shockinitiation behavior 13 can be probed at the atomic scale. The crystal structure of PETN was studied experimentally via xray diffraction by Booth and Llewellyn 11 whose data was later revised by Trotter 12 A tetragonal lattice with spacegroup symmetry cP12 4 was determined, with two C(CH 2 ONO 2 ) 4 molecules (58 atoms) per unit cell 11 Importantly, it was noted that we ak van der Waals forces play an important role in the intermolecular interactions 11 We denote the PETNI polymorph as PETN, but a second polymorph, PETNII, has been reported at high temperature 13 and a third polymorph, PETNIII, has recently been reported near 6 GPa 14 In addition, Olinger et al. 15 performed xray diffraction experiments of linear and volume compression of PETN up to 10 GPa at room temperature. Several theoretical studies of energetic materials including PETN may be found in the literature 1628 Of particular importance for comparison with this work are the firstprinciples studies by Gan et al. 18 Byrd and Rice 23 and Brand 19 Gan et al. 18 performed DFT calculations using the PerdewBurkeErnzerhof (PBE) functiona l and a Gaussian basis set to simulate the hydrostatic compression of PETN. Byrd and Rice 23 have performed a comprehensive study of the accuracy of DFT in applica tion to several EMs. They used DFT with a planewave basis set using the PerdewWang (PW91), PBE, and local density approximation (LDA) flavors of density functional theory at vari ous cutoff energies. Also, Brand 19 performed HartreeFock 39 PAGE 51 calculations of hydrostatic compressi on with a Gaussian basis set. Most of the calculations deal with equilibrium properties or hydrostatic compression. Therefore, the goal of this work is to ex tend the isotropic EOS for PETN to study its anisotropic properties upon hi gh compression. The uniaxial co mpression studies, motivated by the highly anisotropic shocksensitivity be havior of the crystal, were performed to investigate relative changes of physical pr operties upon uniaxial compression along seven low index crystallographic direc tions, <100>, <010>, <001>, <1 10>, <101>, <011>, and <111>. We also completed studies of the equilibrium properties and the hydrostatic EOS of PETNI for the purpose of evaluating the accuracy of DFT as applied to EM crystals and for comparison with previous calculations. The shear stresses upon uniaxial compression are of particular interest because they are usually considered to be the driving forces of plastic deformations in crystals and the sources of mechanochemical reactions behind the shock wave front. Therefore, we investigate the behavior of shear stresses in PETN upon uniaxial compression with the specific aim of seeing if there is any correlation with observed anisotropic shock sensitivity. II. COMPUTATIONAL DETAILS Firstprinciples DFT 29 calculations were performed using the Vienna AbInitio Simulation Package 30, 31 (VASP). To obtain accurate result s, tests were performed to choose appropriate the DFT functiona l, pseudopotential, kpoint set, and energy cutoff the parameters that control the quality of the calculations. The DFT functionals tested were PW91 32, 33 and PBE 34, 35 of the DFT Generalized Gr adient Approximation (GGA). VASP employs either ultrasoft pseudopotentials 36, 37 (USP) or the projectoraugmented wave 38, 39 (PAW) method to eliminate explicit tr eatment of core el ectrons during selfconsistent electronic structure calculations. Therefore, both PW91 and PBE functionals were tested with PAW and USP potentials. For eac h, the experimental structure of PETN was relaxed with an energy cutoff of 1,000 eV and a 2x2x2 MonkhorstPack 40 (MP) kpoint grid. We found that the PBE functional and the PAW potential produced better agreement than any other combination of functional and potential. The energy cutoff is a very important parameter that controls the completeness of the planewave basis set. Together with the density of kpoint samp ling of the Brillouin zone, it determines the numerical accuracy of DFT planewave calculations. It is also known that the DFT calculations of molecular crystals, which are characterized by a highly inhomogeneous distribution of electron density, require larger basis sets to sample the very small electron density in the intermolecular space of the crysta l. Therefore, extensive tests were performed to choose the energy cutoff and kpoint density that would provide sufficient accuracy of total energies, atomic forces, and stresses acting on the PETN crystal both at equilibrium and at highly compressed states. First, total energy calcul ations without relaxation were performed on both the experimental structure and a halfcompressed st ructure using an energy cutoff of 1,000 eV to test for a sufficient kpoint set. The halfcom pressed structure was made by scaling the lattice parameters by 0.5 (1/3) and it was used in this test to be su re that, as the Bril louin zone enlarged during compression, the fixed kpoint set would s till provide results with sufficient accuracy. We found that the 2x2x2 MP set yielded convergence results better than 0.0007 eV/atom in energy, 0.02 eV/A in forces and 0.05 GPa in pressure, and was chosen for the equilibrium structure and compression calculations. It should be noted here that the kpoint spacing in this 40 PAGE 52 grid corresponded to an aver age kpoint spacing of 0.08 1 in each direction when the structure was scaled to half of its experimental volume. Energy cutoff values of 400, 500, 700, a nd 1,000 eV were then tested on the experimental structure via total energy calcu lations using the PBE functional, the PAW pseudopotentials, and the 2x2x2 MP grid. The minimu m energy cutoff, and thus basis set size, with sufficient accuracy was desired. From th is comparison, the 700 eV energy cutoff was chosen to be used for all future calculations, which yielded convergence better than 0.4 GPa in pressure, 0.0015 eV/atom in energy, and 0.015 eV/ in forces. It should be noted that 700 eV is almost twice the nominal energy cu toff for the PAW potential used, 400eV. To obtain the theoretical cell parameters for PETN, the experimental structure was relaxed using the parameters above via the quasiNewton stru cture minimization method as implemented in VASP. Unit cell size and shape, as well as atomic coordinates, were allowed to simultaneously relax without a constraint on symmetry. The convergence criterion for electronic steps for all calculations was set to 10 6 eV, whereas ionic steps were stopped when the maximum force on any atom was less than 0.01 eV/. Hydrostatic compression of PETN was simulated by VASP calculations that compressed the unit cell from V/V o = 1.00 to 0.60 in increments of 0.02, where V o is the predicted volume at zero pressure. At each step within the compression, the unit cell shape and atomic coor dinates were allowed to relax under constant volume without symmetry constraints. The ma ximum force on any atom was less than 0.03 eV/ for all compression simulations. Uniaxial compressions were performed along seven low index crystallographic directions: <100>, <010>, <001>, <110>, <1 01>, <011>, and <111>. For each compression direction, the Cartesian coordinate system, used to specify components of the lattice vectors, was rotated to orient the xaxis along the speci fic crystallographic direc tion to be compressed. Then, the xcomponent of each lattice vector was scaled by increments of 0.02 up to and down to 0/1.06 VV 0/0.7 VV 0 The uniaxial expansions and compressions near equilibrium were also used to obtain elastic co nstants of PETN. At each step, only the atomic coordinates were allowed to relax while unit cell was held fixed. To avoid any jumps in physical quantities due to insufficient sampling of the Brillouin zone, the MP kpoint set was generated for each uniaxial compression to ensure an average grid spacing of 0.08 1 in each direction 41 when the structure was scaled to half of its experimental volume along the corresponding compression direction. III. EQUILIBRIUM PROPERTIES AND HYDROSTATIC EOS The structure of PETNI at zero pressure was calculated, and the unit cell parameters were compared with the data refined by Trotter 12 from the experiments of Booth and Llewellyn 11 Table 1. Calculated equilibrium structur e of PETN compared with experiment. Work a () b () c () V ( 3 ) Booth et al.: expt.[Ref. 11] 9.38 9.38 6.70 589.50 This work 9.617 (+2.5%) 9.612 (+2.5%) 6.826 (1.9%) 630.99 (+7.0%) 41 PAGE 53 see Table 1. The percent error of our calculations of equilibrium lattice constants is within 23%. As in previous DFT studies of energetic materials 17, 21 the prediction of greater lattice constants was believed to be due to an inadequate description of van der Waals dispersive interactions in DFT functionals. Predicted lattice angles, similar to experiment 11, 12 were calculated to be 90.0 degrees. Figure. 1. Isothermal hydrostatic EOS of PETNI. The PETNI isothermal equation of state was obtained by perform ing DFT hydrostatic compression simulations at 0 K. Fig. 1. compares the DFT results with the roomtemperature experimental data from the work of Olinger, Halleck, and Cady 15 The range in volume changes shown is from 60 to 100% of the calculated equilibrium volume, 630.99 3 At V/V 0 = 0.6, the calculated pressure is approximately 30 GPa, which is also approximately the detonation pressure of PETN 42 The calculated pressure at a given volume is higher th an the experimental value, which is believed to be due to the inadequate description of we ak van der Waals dispersive interactions by DFT 21 Meanwhile, the discrepancy between the e xperimental and calculated isotherms is reduced as volume decreases, which is probably due to increased overlap of electron densities of neighboring molecules leading to a better description of th e intermolecular interactions by DFT. This trend has been observed in a number of hydrostatic compression studies of energetic molecular crystals using DFT 21 The lattice constants a and c as a function of pressure are compared with experimental data 15 see Fig. 2. As pressure increase s, the experimental value of a approaches the calculated value. However, the calculated value of c is greater than the experimental value for all pressures. This trend has been observed pr eviously in the DFT calculations of Brand 19 Bond lengths of a PETN molecu le within the crystal were plotted as a function of V/V 0 in Fig. 3. As observed in a previous DFT study by Gan, Sewell, and Challacombe 18 the CC bond and the ON bond experience greater changes th an the other bonds. Th e greatest change in bond length observed was for the proximal ON bond, which decreased by 0.086 A. The bulk modulus, dVdPVB 0 and its derivative with respect to pressure 0B were obtained from fitting the PETNI EOS data to the BirchMurnaghan EOS 18 Figure. 2. Lattice parameters of PETN under hydrostatic compression. Figure. 3. Bond lengths of PETN under hydrostatic compression. 42 PAGE 54 ,)1)(4( 4 3 1)( 2 33 2 0 3 5 3 7 0 B BP (1) where = V/V 0 using pressures from about 0.3 to 10.4 GPa (V/V 0 = 1.06 to 0.72). Following the work of Gan, Sewell, and Challacombe 18 we have restricted our fitting to pressures below 10.45 GPa to be consistent with th e experimental work of Olinger 15 Table 2 shows our data as compared with the BirchMurnaghan fitting 18 of the experimental data obtained by Olinger et al. 15 and the bulk modulus obtained via experimental elastic constants 43 13 33 1211 2 13 12113342 2 CCCC CCCC B (2) from Winey and Gupta 44 Not included in Table 2 are the va lues of the bulk modulus from the theoretical works of Gan et al. 18 (B 0 =14.5 GPa, B' 0 =6.7) which were obtained by fitting their EOS data to the BirchMurnaghan EOS (1) and Brand 45 obtained B 0 =9.38 GPa by using elastic constants via equation (2). Table 2. Bulk modulus results. Work B B 0 (GPa) B B 0 This work (BirchMurnaghan) 9.1 8.3 Olinger et al. (BirchMurnaghan) [Ref. 15] 9.4 11.3 Winey and Gupta [Ref. 44] 9.1 The bulk modulus can be found by fitting to other equations, such as the Murnaghan equation and the Hugon iot conservation equa tions used by Gan et al. 18 The values of the bulk modulus found by fitting our data set to these eq uations were within about 10% of the value obtained by fitting to (1), similar to the results obtained by Gan et al. 18 Using the uniaxial compression data discussed below, the energy E of the crystal as a function of strain was fit to a fourthorder polynomial E( ) The elastic constants calculated were then found from the relation 12 2 0 ii iiE V C (3) where i =1, 2, 3. The calculated values for these el astic constants are in good agreement with the experimental data analyzed by Winey and Gupta 44 see Table 3. Table 3. Elastic constants of PETN. Work C 11 (GPa) C 22 (GPa) C 33 (GPa) This work 18.3 18.5 14.2 Winey and Gupta [Ref. 44] 17.22 17.22 12.17 43 PAGE 55 Figure. 4. Stress tensor components Figure. 5. Change in energy per atom upon uniaxial and hydrostatic compressions. Figure. 6. Band gap upon uniaxial and hydrostatic compression. xx (a), yy (b), and zz (c) of PETN under uniaxial compression. The stresses for the hydrostatic compression are included for comparison to show the degree of anisotropy upon uniaxial compression. IV. UNIAXIAL COMPRESSIONS The diagonal elements of the stress tensor, xx yy and zz, as a function of compression ratio from V/V 0 = 0.7 to 1.0 for each uniaxial compressi on are shown in Figs 4a, 4b, and 4c, respectively. The anisotropic character of cons titutive relationships is clearly seen from Fig. 4a, which shows substantial differences in longitudinal stresses between uniaxial compressions along different crystallographic dire ctions as compared with pressure from hydrostatic compressions when x xyyp zz For strains below 15%, the calculations indicate that the crys tal is less compressible in the < 011> and <101> directions and more compressible in the <110> and <001> directions However, for strains greater than 25%, the <001> and <110> directions emerge as the least compressible. Based upon the work of Jindal and Dlott 8 the anisotropic compressibility of P ETN could be an important factor for determining if anisotropic heating upon shock compression plays a significant role in its anisotropic detonation sensitiv ity. However, we do not calculate the Gruneisen parameters necessary for such a prediction 8 For shocks near 4 GPa in the <110> direction, evidence of detonation has been shown in several experiments 3 Therefore, it is likely that the anisotropic constitutive relati onships can not be probed experi mentally when compressed in this direction beyond 4 GPa. However, the <100> direction has been shown to be very insensitive 13, 5 and therefore it is most likely that the experimental dependence of shear stress on strain can be obtained for this direction in the compression interval shown in Fig. 4. The change in energy per atom and th e band gap of PETNI as a function of 44 PAGE 56 compression for each uniaxial compression and the hydrostatic compression are shown in Figs 5 and 6, respectively. Anisotropic behavior is cl early exhibited in both graphs. The greatest change in energy is experienced when the cr ystal is compressed in the <101> or <011> directions, and all uniaxial comp ressions raise the energy of the crystal more than hydrostatic compression because of the constr aints of the unit cell in transv erse directions. Similarly, the greatest reduction in the band gap is ob served for the <101> and <011> uniaxial compressions. The smallest change in band gap is observed for hydrostatic compression, but it should be noted that the most sensitive direc tion, <110>, experiences only a slightly larger reduction in band gap than the hydrostatic compression. V. SHEAR STRESSES The shear stresses could be of paramount im portance in understand ing the shockinduced anisotropy of EMs. The shockloading cond itions are characterized by fast uniaxial compression of the crystal along a specific crystallographic direction corresponding to the direction of the shockwave propagation. The shear stresses are the driving forces of shearinduced plastic deformations and mechanochem ical events behind the shockwave front. Therefore, we calculate the maximum shear stresses for each uniaxial compression, see Figs 7a, and 7b. Under the uniaxial compression, the crystal experiences the maximum shear stresses x y and x z directed at 45 to the compression direction, their values being 2xyxxyy and 2xzxxzz The behavior of the shear stresses for each compression is relatively similar for the xy and xz maximum shear stresses. The yz values for each compression are relatively different and significantly smaller than the xy and xz For compressions greater than V/V 0 = 0.76, the more sensitive directions, <110> and <001>, show greater magnitude of shear stresses xy and xz than the insensitive direc tions, <100> and <101>. Based upon the greater shear stresses displayed for compression in the more sensitive directions, we have previously suggested a correlation between sensitivity and shear stress behavior under uniaxial compression 4648 Here, we also examine the shear stress maxima, xy and xz as a function of longitudinal stress, xx to compare with experiment and examin e the possibility of a correlation between sensitivity and shear stresses in PETN, see Figs. 8a and 8b. Note that each curve terminates at V/V 0 = 0.70. The calculated longitudinal stress, xx required to compress PETN to the point Figure. 7a. Shear stress maximum xy (a) and xz. (b). 45 PAGE 57 Figure 8a. Maximum shear stress xy (a) and (b) as a function of longitudinal stress xz xx where <110> and <001> compressions exhibit both shear stresses greater than the sensitive <100> and <101> compressions, varies from ~4 GPa for <100> direction to ~12 GPa for <101> direction. As stated above, <110> co mpressions have experienced transition to detonation in Dicks experiment upon compre ssion by longitudinal stresses as low as 4.2 GPa 3 Further, compressions of PETN in the < 001> direction have produced detonation at longitudinal stresses just above 12 GPa 3 The values of xy and xz for the <001> compression do emerge as the greatest values for the compressions studied beyond 12 GPa, while below this value both shear stress maxima for <101> compression remain the gr eatest. However, in Dicks experiment shock in <101> does not lead to detonation for the input stresses up to ~15 GPa, though it exhibits an in termediate velocity transition (abrupt in crease of the shock velocity and stress) at 8.8 GPa. Thus, while there are so me indications of a correlation between maximum shear stresses and shock sensitivity for the <100>, <110>, and <001> cases,, additional analysis needs to be performed. In particular, as suggested in the work of J. Dick 1 it might be helpful to consider the projection of shear stress maxima on available slip systems under uniaxial compression, and this will be the subject of our future work. VI. CONCLUSIONS Firstprinciples simulations of hydrostatic compression and uniaxial compressions in the <100>, <010>, <001>, <110>, <101>, <011>, and <111> directions of PETN have been performed using DFT with the PBE functional and the PAW potential. The calculated lattice parameters were overestimated by about 3%, which is believed to be due to a poor description of van der Waals dispersive in teractions in the functionals of DFT. The calculated bulk modulus and elastic constants are in reasonable agreement with experiment. With the application of hydrostatic pressu re, the calculated un it cell volume and lat tice parameters of PETN show increasing agreement with experimental data. The energy per atom, band gap, and stress es were calculated for each uniaxial compression and compared with the calculated results hydrostatic compression. Our results show a substantial anisotropy of various pr operties of PETNI upon uni axial compressions. The maximum shear stresses were calculate d and examined for a correlation with the anisotropy in shockinitiation sensitivity. The calculated maximum shear stresses, xy and xz upon compression in the sensitive <110> and < 001> directions are greater than for the insensitive <100> direction at input longitudina l stresses above 4 GPa, but no clear correlation 46 PAGE 58 between maximum shear stresses and sensitivity is observed for the <101> direction up to 12 GPa. In the future, we will perform a more deta iled analysis of the sh ear stress behavior under uniaxial compressions by consid ering slip systems of the PETN crystal available upon compression in specific crystallographic direct ions, as suggested by the work of J. Dick 1 to establish more certain correlations with anisotropic shock sensitivity of PETN. ACKNOWLEDGEMENTS The work at USF was supported by the Office of Naval Research (ONR) through the Naval Research Laboratory (NRL) (grants N00173061G022 and N00173082C002) and partly by the Army Research Office through the MultiUniversity Research Initiative on Insensitive Munitions (grant W901NF0510266). The work at Caltech was suppor ted by the Office of Naval Research (ONR) (grant N000140510778) and the Army Research Office through the MultiUniversity Research In itiative on Insensitive Munitio ns (grant W911NF0510345). The work at NRL was also supported by ONR both directly and through NRL. The computations were performed using NSF Teragrid computational facilities (grant TGDMR070018N). REFERENCES 1 J. J. Dick, Applied Physics Letters 44 859 (1984). 2 J. J. Dick, Journal of Applied Physics 81, 601 (1997). 3 J. J. Dick, R. N. Mulford, W. J. Spencer, D. R. Pettit, E. Garcia, and D. C. Shaw, Journal of Applied Physics 70, 3572 (1991). 4 P. M. Halleck and J. Wackerle, Journal of Applied Physics 47 976 (1976). 5 C. S. Yoo, N. C. Holmes, P. C. Souers, C. J. Wu, F. H. Ree, and J. J. Dick, Journal of Applied Physics 88, 70 (2000). 6 L. Soulard, L'Universite de HauteAlsacc, 1990. 7 L. Soulard and F. Bauer, in Shock Compression of Condensed Matter 1989 edited by S. C. Schmidt, J. N. Johnson and L. W. Davison (Elsevier, Amsterdam, 1990), p. 817. 8 V. K. Jindal and D. D. Dlott, Journal of Applied Physics 83, 5203 (1998). 9 P. W. Cooper and S. R. Kurowski, Introduction to the Technology of Explosives (VCH Publishers, Inc., New York, 1996). 10 J. Akhavan, The Chemistry of Explosives (Royal Society of Chemistry Cambridge, 2004). 11 A. D. Booth and F. J. Llewellyn, Jour nal of the Chemical Society, 837 (1947). 12 J. Trotter, Acta Crystallographica 16, 698 (1963). 13 H. H. Cady and A. C. Larson, Acta Crystallographica Section BStructural Science 31 1864 (1975). 14 O. Tschauner, B. Kiefer, Y. Lee, M. Pravica, M. Nicol, and E. Kim, Journal of Chemical Physics 127, 094502 (2007). 15 B. Olinger, P. M. Halleck, and H. H. Cady, Journal of Chemical Physics 62, 4480 (1975). 16 M. M. Kuklja, F. J. Zerilli, and S. M. Peiris, Journal of Chemical Physics 118, 11073 (2003). 17 E. F. C. Byrd, G. E. Scuseria, and C. F. Chabalowski, Journal of Physical Chemistry B 47 PAGE 59 108, 13100 (2004). 18 C. K. Gan, T. D. Sewell, and M. Challacomb, Physical Review B 69, 035116 (2004). 19 H. V. Brand, Journal of Physical Chemistry B 109, 13668 (2005). 20 M. M. Kuklja, S. N. Rashkeev, and F. J. Zerilli, Applied Physics Letters 89, 071904 (2006). 21 H. Liu, J. J. Zhao, D. Q. Wei, and Z. Z. Gong, Journal of Chemical Physics 124 124501 (2006). 22 F. J. Zerilli and M. M. Kuklja, Journal of Physical Chemistry A 110, 5173 (2006). 23 E. F. C. Byrd and B. M. Rice, Journal of Physical Chemistry C 111, 2787 (2007). 24 F. J. Zerilli, J. P. Hooper, and M. M. Kuklja, Journal of Chemical Physics 126, 114701 (2007). 25 F. J. Zerilli and M. M. Kuklja, Journal of Physical Chemistry A 111, 1721 (2007). 26 W. F. Perger, S. Vutukuri, Z. A. Dreger, Y. M. Gupta, and K. Flurchick, Chemical Physics Letters 422, 397 (2006). 27 W. F. Perger, J. J. Zhao, J. M. Winey, and Y. M. Gupta, Chemical Physics Letters 428 394 (2006). 28 W. F. Perger, R. Pandey, M. A. Blanco, a nd J. J. Zhao, Chemical Physics Letters 388 175 (2004). 29 P. Hohenberg and W. Kohn, Physical Review B 136, 864 (1964). 30 G. Kresse and J. Furthmuller, Physical Review B 54, 11169 (1996). 31 G. Kresse and J. Furthmuller, Computational Materials Science 6, 15 (1996). 32 J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh, and C. Fiolhais, Physical Review B 46, 6671 (1992). 33 J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh, and C. Fiolhais, Physical Review B 48, 4978 (1993). 34 J. P. Perdew, K. Burke, and M. Ernzerhof, Physical Review Letters 77, 3865 (1996). 35 J. P. Perdew, K. Burke, and M. Ernzerhof, Physical Review Letters 78, 1396 (1997). 36 D. Vanderbilt, Physical Review B 41, 7892 (1990). 37 G. Kresse and J. Hafner, Jour nal of PhysicsCondensed Matter 6, 8245 (1994). 38 P. E. Blochl, Physical Review B 50, 17953 (1994). 39 G. Kresse and D. Joube rt, Physical Review B 59, 1758 (1999). 40 H. J. Monkhorst and J. D. Pack, Physical Review B 13, 5188 (1976). 41 Materials Studio. 42 L. G. Green and E. L. Lee, Proceedings of the 13th International Detonation Symposium (2006). 43 C. E. Morris, in Sixth Symposium (International) on Detonation edited by J. Short (Office of Naval Research, Arlington, VA, 1976), p. 396. 44 J. M. Winey and Y. M. Gupta, Journal of Applied Physics 90 1669 (2001). 45 H. V. Brand, Chemical Physics Letters 418, 428 (2006). 46 S. V. Zybin, L. Zhang, H. Kim, A. C. T. v. Duin, and W. A. Goddard, in Proceedings of the 13th International Detonation Symposium edited by S. M. Peiris, Norfolk, VA, 2006), p. 848. 47 M. Conroy, I. I. Oleynik, S. V. Zybin, and C. T. White, in Proceedings of the 13th International Det onation Symposium edited by S. M. Peiris, Norfolk, VA, 2006), p. 1191. 48 M. Conroy, I. I. Oleynik, S. V. Zybin, and C. T. White, Shock Compression of Condensed Matter in press. 48 PAGE 60 REFERENCES 1 J. J. Dick, Applied Physics Letters 44 (1984). 2 R. M. Martin, Electronic Structure: Basic Theory and Practical Methods (Cambridge University Press, Cambridge, 2004). 3 M. C. Payne, M. P. Teter, D. C. Allan, T. A. Arias, and J. D. Joannopoulos, Reviews of Modern Physics 64 (1992). 4 P. Hohenberg and W. Kohn, Physical Review B 136 (1964). 5 W. Kohn and L. J. Sham, Physical Review 140 (1965). 6 P. Gianozzi, in Density Functional Theory for El ectronic Structure Calculations 2005). 7 J. P. Perdew and A. Zunger, Physical Review B 23 (1981). 8 J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh, and C. Fiolhais, Physical Review B 46 (1992). 9 J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh, and C. Fiolhais, Physical Review B 48 (1993). 10 J. P. Perdew, K. Burke, and M. Ernzerhof, Physical Review Letters 78 (1997). 11 J. P. Perdew, K. Burke, and Y. Wang, Physical Review B 54 (1996). 12 J. C. Phillips, Physical Review 1 (1958). 13 H. J. Monkhorst and J. D. Pack, Physical Review B 13 (1976). 14 P. W. Cooper and S. R. Kurowski, Introduction to the Technology of Explosives (VCH Publishers, Inc., New York, 1996). 15 J. Akhavan, The Chemistry of Explosives (The Royal Society of Chemistry, 2004). 16 J. J. Dick, Journal of Applied Physics 81 (1997). 17 J. J. Dick, R. N. Mulford, W. J. Spencer, D. R. Pettit, E. Garcia, and D. C. Shaw, Journal of Applied Physics 70 (1991). 18 C. S. Yoo, N. C. Holmes, P. C. Souers, C. J. Wu, F. H. Ree, and J. J. Dick, Journal of Applied Physics 88 (2000). 19 M. D. Segall, P. J. D. Lindan, M. J. Probe rt, C. J. Pickard, P. J. Hasnip, S. J. Clark, and M. C. Payne, Journa l of PhysicsCondensed Matter 14 (2002). 20 G. Kresse and J. Furthmuller, Physical Review B 54 (1996). 21 G. Kresse and J. Furthmuller, Computational Materials Science 6 (1996). 22 G. Kresse and D. Joube rt, Physical Review B 59 (1999). 23 P. E. Blochl, Physical Review B 50 (1994). 49 PAGE 61 24 I. I. Oleynik, M. Conroy, S. V. Zybin, L. Zhang, A. C. v. Duin, I. W. A. Goddard, and C. T. White, in 14th American Physical Society Topical Conference on Shock Compression of Condensed Matter edited by M. D. Furnish, M. Elert, T. P. Russell and C. T. White (American Institu te of Physics, Baltimore, MD, 2005), p. 573. 25 M. Conroy, I. I. Oleynik, S. V. Zybin, and C. T. White, in 13th International Detonation Symposium edited by S. Peiris, Norfolk, VA, 2006). 26 M. Conroy, I. I. Oleynik, S. V. Zybin, and C. T. White, in Shock Compression of Condensed Matter in press. 27 M. Conroy, I. I. Oleynik, S. V. Zybin, and C. T. White, in Shock Compression of Condensed Matter in press. 50 PAGE 62 APPENDIX 51 