USF Libraries
USF Digital Collections

First-principles studies of shock-induced phenomena in energetic materials


Material Information

First-principles studies of shock-induced phenomena in energetic materials
Physical Description:
Landerville, Aaron Christopher
University of South Florida
Place of Publication:
Tampa, Fla
Publication Date:


Subjects / Keywords:
Bimolecular collision
Uniaxial compression
Molecular dynamics
Dissertations, Academic -- Physics -- Masters -- USF   ( lcsh )
non-fiction   ( marcgt )


ABSTRACT: An understanding of the atomic-scale features of chemical and physical processes taking place behind the shockwave front will help in addressing some of the major challenges in energetic materials research. The high pressure shockwave environment can be simulated using computational techniques to predict mechanical and chemical properties of a shocked material. Density functional theory calculations were performed to investigate uniaxial compressions of diamond and both hydrostatic and uniaxial compressions of TATB and NEST-1. For diamond, we calculated shear stresses for uniaxial compressions in the , , and directions and discovered the anomalous elastic regime which is responsible for the significant delay of plastic deformation behind a shockwave. For TATB, the hydrostatic equation of state, bulk modulus, and equilibrium structure were calculated using an empirical van der Waals correction.The principal stresses, shear stresses, and energy change per atom calculated for uniaxial compressions in the directions normal to the {001}, {010}, {011}, {100}, {101}, {110}, and {111} planes show highly anisotropic behavior. A similar study was performed for the newly synthesized energetic material NEST-1 in order to predict mechanical properties under uniaxial compression. From the similarities in the calculated principal stresses for each compression direction we conclude that NEST-1 is likely to exhibit relatively isotropic behavior as compared to other energetic materials. Finally, reactive molecular dynamics of shock-induced initiation chemistry in detonating PETN was investigated, using first-principles density functional theory, in order to identify the reaction mechanisms responsible for shock sensitivities in energetic materials.The threshold collision velocity of initiation for each orientation was determined and correlated with available experimental data on shock sensitivity. The production of NO₂ was found to be the dominant reaction pathway in every reactive case. The simulations show that the reactive chemistry of initiation occurs at very short time scales ~10E⁻¹³ s at highly non-equilibrium conditions, and is driven by dynamics rather than temperature.
Thesis (M.S.)--University of South Florida, 2009.
Includes bibliographical references.
System Details:
Mode of access: World Wide Web.
System Details:
System requirements: World Wide Web browser and PDF reader.
Statement of Responsibility:
by Aaron Christopher Landerville.
General Note:
Title from PDF of title page.
General Note:
Document formatted into pages; contains 68 pages.

Record Information

Source Institution:
University of South Florida Library
Holding Location:
University of South Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 002029127
oclc - 436877672
usfldc doi - E14-SFE0002902
usfldc handle - e14.2902
System ID:

This item is only available as the following downloads:

Full Text
xml version 1.0 encoding UTF-8 standalone no
record xmlns http:www.loc.govMARC21slim xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.loc.govstandardsmarcxmlschemaMARC21slim.xsd
leader nam 2200385Ka 4500
controlfield tag 001 002029127
005 20090916154455.0
007 cr mnu|||uuuuu
008 090916s2009 flu s 000 0 eng d
datafield ind1 8 ind2 024
subfield code a E14-SFE0002902
QC21.2 (Online)
1 100
Landerville, Aaron Christopher.
0 245
First-principles studies of shock-induced phenomena in energetic materials
h [electronic resource] /
by Aaron Christopher Landerville.
[Tampa, Fla] :
b University of South Florida,
Title from PDF of title page.
Document formatted into pages; contains 68 pages.
Thesis (M.S.)--University of South Florida, 2009.
Includes bibliographical references.
Text (Electronic thesis) in PDF format.
ABSTRACT: An understanding of the atomic-scale features of chemical and physical processes taking place behind the shockwave front will help in addressing some of the major challenges in energetic materials research. The high pressure shockwave environment can be simulated using computational techniques to predict mechanical and chemical properties of a shocked material. Density functional theory calculations were performed to investigate uniaxial compressions of diamond and both hydrostatic and uniaxial compressions of TATB and NEST-1. For diamond, we calculated shear stresses for uniaxial compressions in the , and directions and discovered the anomalous elastic regime which is responsible for the significant delay of plastic deformation behind a shockwave. For TATB, the hydrostatic equation of state, bulk modulus, and equilibrium structure were calculated using an empirical van der Waals correction.The principal stresses, shear stresses, and energy change per atom calculated for uniaxial compressions in the directions normal to the {001}, {010}, {011}, {100}, {101}, {110}, and {111} planes show highly anisotropic behavior. A similar study was performed for the newly synthesized energetic material NEST-1 in order to predict mechanical properties under uniaxial compression. From the similarities in the calculated principal stresses for each compression direction we conclude that NEST-1 is likely to exhibit relatively isotropic behavior as compared to other energetic materials. Finally, reactive molecular dynamics of shock-induced initiation chemistry in detonating PETN was investigated, using first-principles density functional theory, in order to identify the reaction mechanisms responsible for shock sensitivities in energetic materials.The threshold collision velocity of initiation for each orientation was determined and correlated with available experimental data on shock sensitivity. The production of NO was found to be the dominant reaction pathway in every reactive case. The simulations show that the reactive chemistry of initiation occurs at very short time scales ~10E s at highly non-equilibrium conditions, and is driven by dynamics rather than temperature.
Mode of access: World Wide Web.
System requirements: World Wide Web browser and PDF reader.
Advisor: Ivan I. Oleynik, Ph.D.
Bimolecular collision
Uniaxial compression
Molecular dynamics
Dissertations, Academic
x Physics
t USF Electronic Theses and Dissertations.
4 856


First Principles Stud ies o f Shock Induced Phenomena i n Energetic Materials by Aaron Christopher Landerville A thesis submitted in partial fulfillment of the requirements for the degree of Master of Science Department o f Physics College of Arts and Sciences University of South Florida Major Professor: Ivan I. Oleynik Ph.D. Brian Space, Ph.D. Sagar A. Pandit, Phd. Dale E. Johnson, Ph.D. Date of Approval: April 10, 2009 Keywords: bimolecular collision, uniaxial compression, molecular dynamics, explosives, detonation, non equilibrium chemistry Copyright 2009 Aaron Christopher Landerville


To my father Paul Euclid Landerville III (1945 2005) Thank you so much for giving me strong legs upon which to carry myself.


ACKNOWLEDGMENTS I would like to thank those with whom I have collaborated in the course of this work. The work for diamond was done in collaboration with Ivan I. Oleynik, Sergey V. Zybin, Mark L. Elert, a nd Carter T. White with support from the Office of Naval Research (ONR) both directly and through the Naval Research Laboratory (NRL ) The work for TATB was done i n collaboration with Mikalai M. Budzevich, Ivan I. Oleynik, Michael W. Conroy, You Lin, and Carter T. White with support from ONR through NRL and partly by the Army Resear c h O ffic e (ARO) through the Multi U niversity Resear c h Initiative on Insensitive Munit ions (Grant No. W901NF 05 1 0266). The work on Nitrate Ester 1 was done in collaboration with Michael W. Conroy, Ivan I. Oleynik, and Carter T. White w ith support from ONR through NRL The work on PETN was done in collaboration with Ivan I. Oleynik and Car ter T. White with support from ONR th rough NRL. All support from ONR and NRL was given by grants ( Grant No. N00173 06 1 G022 and No. N00173 08 2 C002). Calculations were performed using the NSF Te ragrid computational facilities (Grant No. TG DMR070018N). I would like to give special thanks to Mikalai M. Bubzevich for his extensive work on TATB, and Michael W. Conroy for his extensive work on Nitrate Ester 1. Finally, I offer my sincerest thanks and respect to Ivan I. Oleynik for the great lengths to whic h he has assured my success.


i TABLE OF CONTENTS List of Tables iii List of Figures v Abstract vii i Chapter 1 Introduction ................................ ................................ ................................ ..... 1 1.1 Detonation ................................ ................................ ................................ .......... 1 1.2 Hot Spots ................................ ................................ ................................ ............ 3 1.3 Anisotropic Shock Sensitivity ................................ ................................ ............ 4 1.4 First Principles Shock Compression ................................ ................................ .. 5 1.5 Chemical Initi ation ................................ ................................ ............................. 6 1.6 Outline of Fol lowing Sections ................................ ................................ ........... 8 Chapter 2 Uniaxial Compression s of Diamond ................................ ............................. 10 2.1 Computational Details ................................ ................................ ..................... 13 2.2 Uniaxial Compressions ................................ ................................ .................... 14 2 .3 Anomalous Elastic Response ................................ ................................ ........... 17


ii 2.4 Conclusions ................................ ................................ ................................ ...... 18 Chapter 3 Hydrostatic and Uniaxial Compressio n s of TATB ................................ ....... 19 3.1 Computation al Details ................................ ................................ ..................... 22 3.2 Equilibrium Properties and Hydrostatic EOS ................................ .................. 23 3.3 Uniaxial Compressions ................................ ................................ .................... 26 3.4 Conclusions ................................ ................................ ................................ ...... 29 Chapter 4 Hydrostatic and Uniaxial Compression s of Nitrate Ester ............................. 31 4.1 Computational Det ails ................................ ................................ ..................... 32 4.2 Equilibrium P roperties and Hydrostatic E OS ................................ .................. 33 4.3 Uniaxia l Compressions ................................ ................................ .................... 35 4.4 Conclusions ................................ ................................ ................................ ...... 38 Chapter 5 Hypervelocity Collisions in Detonating PETN ................................ ............. 39 5.1 Computational Details ................................ ................................ ..................... 40 5.2 Threshold Velocities t o Initiate Chemical Reactions ................................ ....... 43 5.3 Product Formation ................................ ................................ ........................... 46 5.4 Reaction Pathways ................................ ................................ ........................... 48 5.5 Reaction Dynam ics ................................ ................................ .......................... 50


iii 5.6 Non Equilibrium Chemistry ................................ ................................ ............ 56 5.7 Conclusions ................................ ................................ ................................ ...... 58 Chapter 6 Summary ................................ ................................ ................................ ....... 59 6.1 Shock Compression of Diamond ................................ ................................ ..... 59 6.2 Shock Compression of TATB ................................ ................................ .......... 59 6.3 Shock Compression of Nitrate Ester ................................ ................................ 60 6.4 Hypervelocity Collisions in Detonatin g PETN ................................ ................ 61 List of References ................................ ................................ ................................ .............. 63 About the Author ................................ ................................ ................................ ... End Page


iv LIST OF TABLES Table 1 Lattice parameters and elastic constants for diamond ............................... 1 4 Table 2 Equilibrium properties of TATB unit cell ................................ ................. 2 4 Table 3 Bulk modulus and first pressure derivative for TATB .............................. 26 T able 4 Equilibrium properties of NEST 1 unit cell ................................ ............... 33 Table 5 Pressure dependence of lattice co nstants and volume for NEST 1 ........... 35 Table 6 Bulk modulus and first p ressure derivative for NEST 1 ........................... 35 Table 7 Summary of collision data for PETN threshold cases ............................... 45


v LIST OF FIGURES Figu re 1 Snapshot of ela stic shock response in diamond ................................ ......... 11 Fig ure 2 Unit cells used for diamond compressions ................................ ................ 12 Fig ure 3 Uniaxial data plotted as a function of compression ratio ........................... 1 5 Figu re 4 TATB unit cell at ambient conditions ................................ ........................ 19 Figure 5 Equation of state for TATB ................................ ................................ ....... 24 Figure 6 Lattice parameters under hyd rostatic compression for TATB ................... 25 Figure 7 Principal stresses for TATB uniaxial compression s ................................ .. 27 Figure 8 Compressed TATB latti ce depicting nitrate bending ................................ 27 Figure 9 Change in energy v s. c ompression ratio for TATB ................................ ... 28 Figure 10 Band gap v s. compression ratio for TATB ................................ ................ 28 Figure 11 Shear stresses fo r TATB uniaxial compressions ................................ ...... 29 Figure 12 Equation of state for NEST 1 ................................ ................................ .... 34 Figure 13 Lattice parameters under hydrostatic compression for NEST 1 ................ 34 Figure 14 Principal stresses for NEST 1 uniaxial compressions ............................... 36


vi Figure 15 Change in energy vs. compression ratio for NEST 1 ................................ 36 Figure 16 Simplified depiction of PETN lattice ................................ ......................... 42 Figure 17 Collision geometry of two PETN molecules ................................ ............. 43 Figure 18 Threshold vel ocities vs. collision orientation for PETN ............................ 44 Figure 19 Number of chemical fragments vs. collision velocity for PETN ............... 46 Figure 20 Bond dissociation energies of the PETN molecular complex ................... 48 Figure 21 Reaction dynamics for PETN {001} threshold case ................................ .. 51 Figur e 22 Reaction dynamics for PETN {100} threshold case ................................ .. 53 Figure 23 Reaction dynamics for PETN {110} threshold case ................................ .. 55


vii Fir st Principles Studies of Shock Induced Phenomena in Energetic Materials Aaron Christopher Landerville ABSTRACT An understanding of th e atomic scale features of chemical and physical processes taking place behind the shockwave front will help in addressing some of the major challenges in energetic materials research. The high pressure shockwave environment can be simulated using computat ional techniques to predict mechanical and chemical properties of a shocked material. Density functional theory calcula tions were performed to investigate uniaxial compressions of diamond and both hydrostatic and uniaxial compressions of TATB and NEST 1 For diamond we calculated shear stresses for uniaxial compressions in the <100>, <110>, and <111> directions and discovered the anomalous elastic regime which is responsible for the significant delay of plastic deformation behind a shockwave For TATB, t he hydrostatic equation of state, bulk modulus, and equilibrium structure were calculated using an empirical van der Waals correction The principal stresses, shear stresses, and energy change per atom calculated for uniaxial compressions in the directions normal to the {001}, {010}, {011}, {100}, {101}, {110}, and {111} planes show highly anisotropic behavior. A similar study was performed for the newly synthesized energetic material NEST 1 in order to predict mechanical properties under uniaxial compressi on. From the similarities in the calculated


viii principal stresses for each compression direction we conclude that NE ST 1 is likely to exhibit relatively isotropic behavior as compared to other energetic materials. Finally, r eactive molecular dynamics of sho ck induced init iation chemistry in detonating PETN was investigated using first principles density functional theory in order to identify the reaction mechanisms responsible for shock sensitivit ies in energetic materials The threshold collision velocit y of initiation for each orientation was determined and correlated with available experimental d ata on shock sensitivity The production of NO 2 was found to be the domi nant reaction pathway in every reactive case The simulations show that the reactive c hemistry of initiation occurs at very short time scale s ~10 13 s at highly non equilibrium conditions and is driven by dynamics rather than temperature.


1 1.0 INTRODUCTION An energetic material (EM) is any material with a large amount of stored chemical energy which can be released with relative eas e. One particular class of EM, known as insensitive munitions, encompasses the many explosives used in both indust rial and military applications. These materials are predominantly composed of molecular crystal s which are insensitive to shock, making them relatively safe to use. Such materials do not release the bulk of their stored energy by thermal decomposition otherwise known as deflagration ; rather they undergo a process called detonation which is induced by delivering a strong shock to the material causing the propagation of a shockwave The major goal of this work is to elucidate the complex physical and chemical processes that occur immediately behind the shockwave front during the initiation of detonati on in order to understand the mechanisms responsible for the shock sensitivity of a material 1.1 DETONATION Detonation Vielle, Mallard, Le Chatelier, and Berthelot began to study flame pro pagation. 1 They noticed that when a flammable gas in a hollow tube was ignited at one end, the flame propagated to the other end of the tube at a rate of only a few meters per second (deflagration). There were exceptions, however, where under certain circu mstance the


2 flame would suddenly begin propagating at rates in excess of 2 km/s This newly discovered form of chemical decomposition was then called detonation. The two completely different types of decomposition were observed again and again in both gase ous and solid explosives, but the reason for the existence of detonation remained a mystery until two physicists named Chapman and Jouguet independently arrived at the same explanation in 1899 and 1905 respectively. 1 They concluded that in the case of deto nation a shockwave propagates through a material and causes ins tantaneous chemical reactions at the high pressure point, called the C J point, which resides at the shockwave front. Subsequently, the Chapman Jouguet process for de tonation beca me the founda tion upon which later models would be built. It is for certain that chemical reactions do not just occur instantaneously; some mechanism over some time scale must be involved. In 1943, a refinement to the Chapman Jouguet model was conceived by Zeldovich von Neumann, and Doering. This model, termed the ZND model of detonation, corrected the discontinuity of the unreacted material and detonation products in the C J model by introducing a zone of reaction. In the ZND model, a high pressure point, called t he von Neumann spike, passes over the material and immediately compresses it. Then, the highly compressed material begins to under go chemical reaction s which cease once the lower pressure C J point is reached. Finally, the reacted material, which is in a gaseous state, expands behind the C J point in a zone of rarefaction. The energy released in the reactive zone between the von Neumann spike and the C J point drives the shockwave forward to produce a steady state detonation process. While the ZND model g ives a more realistic view of an actual


3 detonation process, it still does not account for the complex chemical mechanisms which occur in the reaction zone. 1.2 HOT SPOTS There is much more to a real detonation process than for which the ZND model can accou nt. In particular, it was found that the formation of defects tends to enhance the shock sensitivity of a material. When a shockwave passes over s mall imperfections in the crystalline environment, such as small voids or crack s, zones of high temperature ca lled hot spots can form giving rise to chemical initiation These hot spots can form either by void collapse or by friction over slip planes within the crystal; in either case a hotspot will grow if the heat being released by chemical reac tion is greater t han the rate at which heat can dissipate into the surrounding crystalline environment 2 Also, the size of a defect contributes greatly to initiation; if the defect is below some critical size, then a steady state detonation p 3 A more rec ent study computationally modeled two dimensional v oids and cracks using molecular dynamics (MD) and found that the amount of ejected material and the temperature resulting from collisions greatly depended upon the size of the defect 4 Their findings sugge sted that the temperature resulting from a passing shockwave increases as a function of gap distance, but planes off once some critical gap distance is reached An interesting counterpoint to the critical defect size is the findings of Bowden and Singh tha t single crystals have a critical minimum in size for detonation to occur that is similar to the critical size of defects for a given material 2 If an EM crystal is hit with an anvil, or heated sufficiently it will detonate provided it is above some critic al size. If it is below this critical size, or is of ample size but heated or


4 impacted insufficiently, it will merely deflagrate and crack along preferred crystallographic planes. This discovery has a couple of importan t implications for detonation. First if the crystal is without defects, then defects can still be formed upon application of shock and heating, t hus still allowing for enhanced sensitivity. Second, the fact that insufficient heating or loading of a crystal will lead to cracks along preferre d crystalline planes is the first glimpse into the importance of anisotropic sensitivities in EM crystals. In fact, a little over thirty years after the work of Bowden and Singh, Jerry Dick observed the anisotropic shock response of PETN at Los Alamos, th us beginning a new chapter in the investigation of detonation. 1.3 ANISOTROPIC SHOCK SENSITIVITY While studying the shock initiation of PETN single crystals, Jerry Dick found that the amount of shock needed to initiate the detonation process depended upo n the crystallographic direction in which the shock was delivered 5 The subsequent work of Dick et al found that PETN was most sensitive in the directions normal to the {001} and {110} planes, and least sensitive for {100} and {110}. 5 7 To explain the obs erved anisotropy, Dick proposed a model of steric hindrance. In this mod el 6,8 he proposed that upon the application of shock, slip systems form along crystallographic planes oriented about the direction of the ma ximum resolved shear stresses. These slip planes are hindered for sensitive directions, and relatively unhindered for insensitive directions The hindrance then gives rise to intermolecular collisions along the slip plane thus creating hot spots which lead to initiation and subsequent detonation. A computational st udy was performed by Christine Wu et al. to test this. 9 In this study, they deformed PETN


5 molecules as would be expected upon compressions in the <001>, <100>, and <110> crystallographic directions They then collided the deformed molecul es along the directions of slip using reactive MD and looked for chemical react ions. They found that for collisions corresponding to shock compressions in the <001> and <110> directions a lower collision velocity was required for chemical initiation leadin g to the conclusion that those directions have enhanced chemical sensitivity. Likewise, they found that collisions corresponding to shock compression in the <100> direction had very low chemical sensitivity By mirroring the experimental results of Dick e t al t he computational experiment of Wu et al. lent further credence to the concept of steric hind rance. 1.4 FIRST PRINCIPLES SHOCK COMPRESSION The steric hindr ance model of detonation has prompted other studies regarding the mechanical properties of EM s. In a recent study done by Mike Conroy et al the principal and shear stresses were inv estigated for PETN 10 Using densi ty functional theory, they performed uniaxial compressions of the lattice vectors of the bulk PETN crystal and allowed the positions o f the atoms to relax. From the relaxed structure at each compression step, they calculated the princ ipal and shear stresses, the change in ener gy per atom, and the band gap They found that the principal and shear stresses as a function of compression we re much greater in the sensitive <001> and <110> directions than in the insensitive <100> direction. With the assumption that the principal and shear stresses can be correlated with shock sensitivity, it seemed that the work of Conroy et al. was yet anothe Though they could not draw a direct correlation between the shear stresses and


6 sensitivities, the work of Conroy et al hinted at a relationship between the physical properties and the sensitivities of EM bulk crystals 1. 5 CHEMICAL INITIATION While it is clear that the mechanical properties of an EM are important for initiation, the importance of the chemical properties should not overlooked. After all, the chemical energy stored in an EM is within the chemical bonds which are broken and reformed during detonation. Since the conceptualization of the Chapman Jouguet process for detonation, physical scientists have been looking for ways to describe the internal chemical mechanisms t hat give rise to init iation at the shockwave front. This prompted the implementation of classical transition state theory first developed by Henry Eyring in 1935. 11 Classical transition state theory provides a good, accurate framework for reaction rates an d pathways using Equation 1 w here f is the mode frequency, E a is the activation energy of the chemical bond associated with the mode, T is the temperature of the system, and k B constant. One of the main assumptions behind classical transition state theory is that there is quasi equilibrium between the activated complexes at the transition state and the reactants. This, along with Equation 1, led to the developmen t of the Eyring equation, which can be used to predict the reaction rate of a chemical process that is in quasi equilibrium. However, the impl ementation of classical transition state theory to detonation processes failed as it became clear that detonating materials react at a slower rate than the predicted rate given the temperat ure of the process 12 This prompted the


7 development of newer models t hat could describe the non equilibrium chemistry of a steady state detonation process es Such efforts were sp earheaded by Eyring hi mself as seen in the development of starvation kinetics 12,13 which proposes that the slowed reaction rate is due to the notion that a breaking bond must draw its energy from a vibrational reservoir that is out of equilibrium with the reacting molecule. 12 Such an environment could presumably be seen in the high pressure zone immediately behind the shockwave front i n which molecules are compressed over extremely small timescales and would involve the transfer of energy from intramolecu lar to intermolecular degrees of freedom T he work of Coffey and Toton, 14 expanded upon the idea of starvation kinetics by showing that multi phonon excitation of the intermolecular vibrational degrees of freedom within a reacting molecule should be possib le. Zerilli and Toton, 15 expanded upon this work by allowing phonon frequencies to vary as a function of pressure. Drawing from Eyring, Coffey, Zerilli, and Toton, Walker offered a similar but alternate view of how energy was passed from intramolecular to intermolecular modes, as well as from intermolecular to other intermolecular modes which was termed physical kinetics 16,17 He proposed that the limiting factor of energy mediation was not just the frequency of intermolecular vibrations, but also the veloc ity of the vibrating atoms themselves; noting that the speed of the passing shockwave was experimentally determined to be of the same order of magnitude as the velocity of the vibrating atoms. 16 In his model, Walker showed that the similarities in the shoc kwave speed and passing time to atomic speeds and vibration periods within a reacting molecule could lead to mechanical, not thermal, bond scission. Finally, Dana Dlott proposed a model of vibrational cooling and multi phonon up


8 pumping to explain the tran sfer of vibrational energy in energetic solids. 18 For this transfer of energy, Dlott introduced the concept of doorway modes which are low frequency intermolecular vibrations which can be excited by phonons from the surrounding environment Subsequently, D lott expanded upon the multi phonon up pumping model in a series of works 18 21 to produce perhaps one of the most successful theoretical treatments of shockwave phenomena in energetic solids to date. But w hile the current models of detonation do well in pr edicting the reaction rates of the steady state detonation process, they fall short of explaining the dynamics and reaction mechanisms which give rise to initiation. 1.6 OUTLINE OF FOLLOWING SECTIONS In the following secti ons, we present four studies that are mutually related to shockwave phenomena. The first three studies involve compressing the material of interest uniaxially to investigate mechanical stress and binding energy on the atomic scale in the highly compressed environment behind the shockwave front. The first study involving diamond tests the accuracy of the ato mistic REBO potential at high compr ession using DFT. Also the validity of the anomalous elastic response seen in shockwave simulations of diamond using MD with the REBO potential is te sted. The next study i nvolves hydrostatic and uniaxial compression s of the very important EM TATB to determine anisotropic constitut ive relationships and the direction al dependence of shock sensitivities This study is very useful in showing the effective ness of DFT with an empirical van der Waals correction in predicting accurate properties of a material. The third study involves hydrostatic and uniaxial compressions of a newly synthesized EM


9 called NEST 1. Little is known about this new material as very little experimental and theoretical work has been done. We predict mechanical properties such as the hydrostatic EOS and the bulk modulus, and go further by making predictions about the anisotropic behaviors and shock sensitivities of NEST 1. Finally, th e last study presented here is an MD study on bimolecular collisions of the very well known EM PETN using DFT. Assuming that the passing of a shockwave produces hypervelocity molecular collisions in energetic solids, we reduce the crystalline environment t o a system of two colliding molecules, which is a manageable size for highly accurate DFT calculations. This study aims at elucidating the mechanisms that lead to the chemical initiation of detonation. We investigate chemical sensitivities by varying the speed of collision and determining threshold velocities of initiation. We also investigate the dissociation energies of various bonds in the system and trace the trajectories of reactants and products at femto second resolution. Despite the differences in materials and techniques for the aforementioned studies, the underlying theme is simple. We aim to elucidate the complex atomic scale processes that occur in the environment immediately behind the shockwave front using highly accurate DFT calculations. A n understanding of these processes will be extremely useful in addressing the major issues that confront researchers in the energetic materials community, particularly the factors that influence the sensitivities of energetic materials to heat and impact s timuli.


10 2.0 UNIAXIAL COMPRESSION S OF DIAMOND Diamond is a unique material with many exceptional physical properties 23 The diamond crystal is very resistant to shear deformation because of a strong angular dependence in the interatomic interactions. In fact, it is one of a handful of materials where the shear moduli, and are larger than the bulk modulus, which is the highest of all elemental s olids. The unusual hardness of diamond is a direct consequence of its shear stiff ness; a property successfully exploited in the diamond anvil cells used as the primary tool for generating ultra high pressures in the l aboratory 24,25 One of the attractive methods to probe extreme mechanical proper ties of diamond is to sub j ect it to shock compression 26 Experimental results for mechanical properties of shock compressed diamond are only beginning to emerge 27 31 Recently, the hydrocarbon reactive bond order (REBO) potential 32,33 was used to perf o r m large scale MD simulations of shock wave propagation in diamond along the < 100 > < 110 > and < 111 > d irections 34 These simulations aimed at stimulating future experimental e ff orts and to p rovide insight into the orientational depend ence of sho ck induced chem istry 35,36 With increasing shock wave intensity, four diff erent regimes of materials r esponse were observed : a pure elastic regime, shock wave splitting into elastic and plastic waves, an


11 anomalous elastic regime, and an overdriven plastic wave with acti vated carbon chemistry 34,36 The anomalous elastic response in shock compressed diamond is characterized by the absence of plastic deformations, including point and extended defects in a highly compressed diamond lattice as depicted in Fig ure 1. The shock wave intensities corresponding to the anomalous elastic regime lie between those resulting in plastic deformation and those that initiate carbon chemistry in highly compressed diamond. Zybin et al. w ere able to relate this regime to a non monotonic behav ior of shear stress as a function of u niaxial compression 36 However this correspondence was made using a hydrocarbon R EBO potential, which is known to provide a good description of mechanical prop erties o f diamond at equilibrium, but has not yet been tes ted at very high pressures and stresses. Therefore, it is unclear whether the anomalous elastic regime found in the MD s imulations is a real phenomenon or just an artifact of this potential. Figure 1 : Snapshots of the shock front in the anomalous elastic region for a <110> shockwave (top) resulting in a compression ratio of 0.795, and a <111> shoc kwave (bottom) resulting in a compression ratio of 0.715. Atoms are colored according to potential energy (blue corresponds to low, red to high energy).


12 In this work, we give results of a DFT investigation of diamond under uniaxial compression at zero t emperature in the < 100 > < 110 > and < 111 > directions in the strain interval We fi nd that the non monotonic behavior of the shear stress under uniaxial compression originally observed in calculati ons using the REBO potential is also seen in the DFT calculations indicating that this is a real response. The results show, however, that either the hydrocarbon REBO potential will have to be modi fi ed or new potentials introduced to better reproduce the high pressure DFT results. The unit cells used in these studies are shown in Figure 2. Di ff erent strains were imposed by varying the c axis of a particular cell. For each value of c the lateral d imensions of the cell were held fi xed and all the atomic c oordinates were relaxed to have zero net forces on the atoms. The starting atomic coordinates were obtained from the atomic coordinates of the previously relaxed structure with larger c Uniaxial compression at fi xed lateral dimensions corresponds to the c onditions of shock wave experiments at su ff iciently large shock wave intensities. The time scale associated with the initial process of shock wave compression is on the order of picoseconds or less (the thickness of the shock wave front is of the order of several lattice constants and the shock wave velocity is several km/s ). Therefore, the lattice rapidly transforms at the shock wave front to a uniaxially compressed state. Figure 2 : Unit cells used in calculating the effects of uniaxial compression along the <100>, <110>, and <111> directions.


13 This condition of pure uniaxial compression of shocked diamond diff ers from conditi ons resulting from indentation or diamond anvil cell compression where the mechanical loading is comparatively slow so that the stresses in the lateral dimensions are able to relax (at least partially) through lateral expansion of the lattice. Hence, our i nitial conditions diff er from those used in several other DFT calculations of shear d eformations, where the conditions of diamond anvil cell experiments (speci fi c combination of uniaxial and hydrostatic strains) 37,38 or conditions of nano indentation a nd d iamond cleavage experiments (full relief of stresses in lateral dimensions) 39 41 were used. Calculations of pure uniaxial compression of diamond similar to our c alculations were performed by Nielsen 42 However, his study did not reach the high degree of co mpression where the shear stresses exhibit non monotonic behavior. 2.1 COMPUTATIONAL DETAILS Our DFT calculations were performed using the Vienna ab initio simulation package (VASP) 43,44 within a generalized gradient approximation (GGA) for the electronic exchange and correlation (PBE GGA functional was used). 45 The projector augmented wave (PAW) potential and a large plane wave cuto ff energy of 700 eV were used. Recent work has demonstrated that even at high compression the pseudopotential and PAW methods give almost indistinguishable results from those ob tained by all electron methods 46 A possible appearance of metallic phases in the course of uniaxial compression requires dense sampling of the k space Brillouin zone, thus the k point spacing was 0.05 1 The particular values of the energy cuto ff and k space sampling density were chosen to achieve accuracy in the calculated stresses of better than 0 1 GPa


14 and in the calculated energies of better than per atom. The calculated value s of lattice and elastic constants are in good agreement with experiment, as shown in Table 1 Table 1 : DFT and REBO bulk properties (lattice parameter in and elastic constants in MBar) of diamond as compared with experiment. Met hod a C 11 C 12 C 44 B C' DFT 3.530 11.0 1.19 5.86 4.42 4.90 REBO 3.556 10.7 1.00 6.8 4.23 4.85 Experiment 3.567 10.7 1.25 5.77 4.43 4.76 2.2 UNIAXIAL COMPRESSIONS We compare DFT to REBO results for the binding energy per atom as a function of physical s train at the top of Fig ure 3. The DFT and REBO binding energies are in good agreement for strains up to 0.25 in the < 100 > and < 111 > directions, and 0.15 in the < 110 > direction. The REBO results plotted in Fig ure 3 cover a smaller strain range than the DFT results because at very large strains we faced problems due to the fi nite cuto ff radius of interatomic interactions. Discontinuities in the slope of the binding energy for the < 110 > and < 111 > directions are the signatures of polymorphic phase transitions o ccurring at large compression ratios. For example, at the < 111 > compressed diamond crystal structure transforms to the structure with atomic arrangements close to those in a simple hexagonal lattice.


15 Figure 3 : Binding energy (top), longitudinal stress (middle), and shear stresses (bottom) for uniaxial compressions along <100>, <110>, and <111> directions. The longitudinal stress, zz as a function of uniaxial strain is shown in Fig ure 3 (middle). G ood agreement between the DFT and REBO data is now observed in a narrower interval of uniaxial strain: for the < 100 > and < 111 > directions, and for the < 110 > direction, which refl ects the fact that the stresses are fi rst order derivatives of energy per unit volume on strain. The longitudinal stress shows monotonic behavior as a function of strain in the < 100 > direction but exhibits a local maxi mum and minimum in the < 110 > and < 111 > directions. Results from DFT and


16 REBO calculations show qualitatively similar behavior, but the positions of the extrema and th e values of the stresses are diff erent. These local extrema are the primary sources of th e non monotonic behavior of shear stresses for the < 110 > and < 111 > directions shown in Fig ure 3. In contrast, the non monotonic behavior of shear stress in the < 100 > direc tion shown at the bottom of Figure 3 is due to the change of slope of both lateral x x(yy) and longitudinal zz stress pro fi les. Although xx (yy) are monotonic with zz > xx(yy) at small strains, they fi rst approach and then becomes larger than zz in the intermediate r ange of Upon further increase of the strain, the slope of zz beg ins to increase resulting in the next crossover, zz = xx(yy) with zz again larger than xx(yy) at very large compressions. Uniaxial compression in the z direction creates the maximum shear stresse s and th a t are directed at 45 to the direction of the compression. It is the shear stress that drives plastic deformations when exceeding a threshold value. The DFT and REBO shear stress pro fi les are shown in Fig ure 3 (bottom). The shear stresses x z and yz are equal for the < 100 > and < 111 > directions because of the symmetry of the diamond lattice. Therefore, only one of them is shown in the left and right panels of Fig ure 3 (bottom). Both results from the REBO potential and DFT exhibit non monotonic behavior of the shear stresses for all three crystallographic directions. The DFT shear stresses initially increase, reach a maximum, decrease, reach a minimum and increase again at very large compressions. The DFT maximum and subsequent minimum are at compression ra tios 0 65 and 0 425 for the < 100 > direction, 0 75 and 0 675 ( x z ) 0 725 and 0 675 ( yz ) for the < 110 > direction, and 0 675 and 0 525


17 for the < 111 > direction, see Fig ure 3 (bottom). The REBO results show similar behavior but due to problems with the fi n ite cuto ff radius, the shear pro fi les are limited to smaller strains. 2.3 ANOMALOUS ELASTIC RESPONSE Both DFT and the REBO potential predict that the shear stresses approach small values for both the < 110 > and < 111 > directions which were studied in the MD shock simulations of McLaughlin et al 34 This causes important changes in the mechanical properties of shock compressed diamond. N o plastic deformations were observed in the MD simulations of McLaughlin et al. in just this range of uniaxial compressions 3 4 I n contrast, in the regime of smaller compressions, plastic deformations occurred by displacement of the atoms in the direction of the maximum shear stress. This causes stress relief in the compressed lattice. The same picture is observed at higher uniax ial compressions when the shear stress again becomes appreciable. Because the shear stress is the driving force for the lateral movement, slipping of the crystal planes, and the creation of both point and extended defects, there will be some range of uniax ial compressions where these processes are inhibited or even frozen due to small values of shear stress near the minima of the shear pro fi les in Fig ure 3 (bottom). This anomalous elastic state might ultimately transform to a plastic state behind the shock front with this transformation initiated by random thermal movement in the s hock heated crystal. However, this instability was not observed in the MD simulations of McLaughlin et al 34 suggesting that this anomalous elastic state, even if unstable, could persist for long enough times to be experimentally observable.


18 2.4 CONCLUSIONS In spite of quantitative problems with the REBO potential at large compressions, the DFT results indicate that the anomalous elastic regime seen in MD simulations using this potential is a real phenomenon that should be observable in su ffi ciently resolved experiments. These phenomena may also occur in several other materials that exhibit non monotonic shear stress strain relationships. 47 Our results also show that further work is required to develop robust and transferable interatomic potentials for a quantitative description of strong shock waves in condensed matter even for a relatively simple system such as diamond. Moreover, our results provide benchmarks for testing other classical carbon potentials at conditions appropriate for strong shock wave simulations.


19 3.0 HYDROSTATIC AND UNIAXIAL COMPRESSION S OF TATB 1,3,5 triamino 2,4,6 trinitrobenzene (TATB) is a high explosive known for its high insensitivity 48 which has be en observed in several standard tests. 48,49 For example, TATB does not react when subjected to friction tests, spark tests, or drop tests. During a cook off test, TATB undergoes a burning reaction, but not detonatio n Despite its high stability and relati ve safety, TATB is a highly powerful organic explosive which makes it an interesting target for both theoretical and experimental studies. The structural and mechanical properties of energetic materials (EMs) such as TATB are important for developing rob ust models of detonation. TATB is an organic molecular crystal with two molecules per unit cell. The unit cell is triclinic with symmetry and has the following lattice parameters: , and 50 The TATB unit cell is shown in Figure 4 Unlike for several other EMs such as PETN, HMX, and RDX, no polymorphs of TATB have been discovered. Infrared spe ctroscopic studies of TATB at high pressures have shown that there are no phase transitions under hydrostatic compressions up to 10 GPa. 51 Figure 4 : Triclinic unit cell of TATB at ambient conditions.


20 Using density functional theory (DFT) calculations, Manaa reached the insulator metal transition in TATB by both ben ding a nitro group out of the molecular plane by an angle of 55 and by increasing the N O bond by 0.1 52 An optical study performed by Cady and Larson revealed a layered structure in which molecules are arranged in planar sheets almost parallel to the a b plane, with an interplanar distance of 3.3 50 It has been suggested in the work by Wu et al. that the intermolecular hydrogen bonding in TATB contributes to its high insensitivity. 53 Further, the work by Pravica et al. suggested the strengthening o f intermolecular hydrogen bonding with pressure. 51 Like other organic EMs, 10 ,54,55 TATB may also exhibit anisotropic sensitivity to shock. A well kn own example of studies involving the anisotropic sensitivity to shock induced detonation is the work of D ick et al ., who first discovered this property in PETN. 54 Dick et al found h igh sensitivity to shock in the <001> and <110> directions, and low sensitivity to shock in the <100> and <101> directions. A recent DFT study of PETN performed by Conroy et al 10 found greater values of shear stress under simulated uniaxial compression for the more sensitive directions found in the work of Dick, suggesting a possible correlation between shear stress and sensitivity. We employ a method similar to that used by Con roy et al to investigate the possible existence of anisotropic shock response in TATB. Unfortunately, there are no experimental data concerning the uniaxial shock response of TATB. However, first principles DFT methods make possible the study of electro nic and mechanical properties of EMs under various types of mechanical loading such as uniaxial compression. For example, Wu et al performed a DFT study on TATB


21 which indicated that band gap closure occurred near 47% compression when uniaxial strain was applied in the c direction indicating that the lower bound for the metallization pressure was 120 GPa. 53 It was argued in Ref. [53 ] that this result does n ot support the works of Williams 56 and Gilman 57 which suggest that electronic excitation, i.e. meta llization, at high pressure plays a major role in the initiation of chemistry and the propagation of detonation waves in EMs. Pure DFT calculations of organic molecular crystals such as EMs suffer from limitations in accuracy. Byrd et al 58 showed that DFT calculations have the tendency to overestimate the equilibrium volume of the unit cell for EMs due to a deficiency in the description of weak hydrogen and van der Waals interactions. 58,59 For example, pure DFT calculations for TATB overestimate the equ ili brium volume by 14% 60 One approach to correct this problem is via first principles methods, requiring nonlocal approximations to account for dispersion interactions (vdW or weak hydrogen forces) within DFT. 61 64 These methods show great promise for t he future of first principles calculations, but the application of these methods at present is too costly for practical calculations of large systems such as EMs. Another approach employs empirical potentials and adds only a negligible cost to DFT calcula tions. 65 68 This method combines pure DFT calculations with an empirical van der Waals correction, which is defined as the sum over atom atom pair wise potentials. The pair wise potentials are formed by a product of the asymptotic te rm for the interaction at large distances and a damping function that both prevents divergence of the potential at large distances and reduces the influence of the empirical potential at distances where pure DFT provides an accurate interaction


22 description. This work utilizes an empirical approach based on the work of Neumann and Perrin 67 that was recently impleme nted in a DFT study of EMs upon compression. 60 The goal of this work is to use vdW corrected DFT calculations to investig ate structural and mechanical properties, as well as the electronic structure, of TATB under hydrostatic and uniaxial compressions normal to the {100}, {010}, {001}, {110}, {101}, {011}, and {111} crystallographic planes. Our particular focus is on the she ar stresses resulting from uniaxial compressions, which are known to play a key role in the development of plastic deformations. We hypothesize that plastic deformations in EMs might be responsible for mechanically inducing chemical reactions that lead to detonation. Further, the previous work of Conroy et al. 10 suggests a possible relation between shear stress and detonation sensitivity. Using this assumption, our goal is to identify directions which might exhibit greater sensitivity in TATB. 3.1 COMPUTA TIONAL DETAILS We used the Vienna Simulation Package (VASP) 43,44 to perform DFT studies of TATB. Th e PW91 functional 69,70 and projector augmented wave (PAW) potentials 71,72 were used in our calculations because the parameters of the empirical vdW corre ction 67 were determined with these settings. Tests similar to those of Conroy et al. 10 ,73 were performed to determine other calculation parameters that would provide an appropriate level of accuracy. An energy cutoff of 700 eV was used, which yielded the energy per atom to within 2 meV. Monkhorst Pack 74 k point sets were chosen for each manner of compression simulation such that the average k point spacing was at most 0.08 1 at a compression of The conjugate gradient structure mi nimization method was


23 used to perform structure optimizations of the TATB crystal. Using these calculation parameters, we found that the conjugate gradient relaxation method reduced the calculation time as compared to the quasi Newton relaxation algorithm. Self consistent field electronic structure calculations were done with a tolerance of 10 6 eV. Hydrostatic compression calculations were performed within the range from to where is th e calculated equilibri um volume, in increments of 2% At each compression step, the lattice parameters and the atomic coordinates were relaxed at constant volume until the maximum force on any atom was less than 0.03 eV/ Uniaxial compressions were also p erformed along seven low index directions which are defined to be the directions normal to the crystallographic planes given by the Miller indices {001} {010} {011} {100} {101} {110} and {111 }. For each uniaxial compression orientation, the unit cel l was rotated such that the compression direction was aligned with the x axis. Subsequently, the x components of each lattice vector were rescaled at each stage of compression. The lattice parameters were held constant such that only the atomic coordinates were allowed to relax 3.2 EQUILIBRIUM PROPERTIES AND HYDROSTATIC EOS The unit cell parameters of TATB were relaxed in order to find the calculated zero pressure structure. In Table 2 the results of the calculation are compared with the experimental st ructure from the data of Cady and Larson and Stevens et al 50 ,75 The a b and c lattice constants are in good agreement with experimental data, yielding less than 1% error as compared to both experiments.


24 Table 2 : Equilibrium u nit cell properties of TATB with comparison between Study a () b () c () Vol. ( 3 ) Cady and Larson 9.010 9.028 6.812 108.59 91.82 119.97 442.49 Stevens et al. 8.967 9.082 6.625 110.5 93.0 118.9 442.5 This Work 9.0231 9.0415 6.585 106.85 92.02 119.75 435.90 % Er ror to Cady and Larson 0.15% 0.15% 3.33% 1.6% 0.22% 0.1 8% 1.49% % E rror to Stevens et al. 0.63% 0.45% 0.58% 3.3% 1.05% 0.72% 1.49% The isothermal equation of state (EOS) for TATB was calculated by performing hydrostatic compressions. The calculated EOS is compared both with experimental data 75,76 a nd theoretical data 59,77 in Figure 5 Our calculated EOS agrees well with Stevens et al and Olinger et al ., but significant disagreement is shown with the theoretical predictions of Pastine and Bernecker, and Byrd and Rice. In addition to the EOS, the la ttice constants a b and c as a function of pressure were also calculated and compared with experiment in Figure 6 The results are in good agreement with experimental data. 76 There were some uncertainties in the calculation of the lattice vectors in the work by Olinger and Cady. 76 In their study, the lattice constants Figure 5 : Isothermal hydrostatic EOS of TATB in the range of V/V 0 = 1.00 0.60.


25 a b and c were not directly measured, but were determined using some assumptions: 1) the a/b ratio remained constant, and 2) any point on a molecule will remain in the same relative positi on with respect to a corresponding point on a neighboring molecule in a adjacent plane. Our results indicate that these are reasonable assumptions. Under hydrostatic compression, the a/b ratio changes in value by less than 0.1% up to In addition, we found that the molecules on two adjacent planes rotate by a maximum angle of 0.07 with respect to each other under a pressure of 27.5 GPa, which we believe to be negligible. Figure 6 : Lattice parameters of TATB under hydrostatic compressi on. We obtained the bulk modulus and its derivative with respect to pressure by fitting the hydrostatic compression data to the Birch Murnaghan EOS 83 wh ere To be consistent with experimental results, our data was fitted up to 8 GPa. Table 3 exhibits our results as well as values report ed in other works and determined from the reported data. 55,75 77 Equation 2


26 Table 3 : Bulk modulus and its pressure derivative for TATB. All values reported were determined by fit to the Birch Murnaghan EOS up to 8 GPa. Study Type (GPa) This Work Theory DFT vdW 13.8 9.4 Olinger and Cady E xperiment 16.7 5.7 Stevens et al. Experiment 13.6 12.4 Byrd and Rice Theory Pure DFT 2.7 16.8 Pastine and Bernecker Theory 6.2 47.7 Reasonable agreement was found for the bulk modulus between our result and the experimental results of Olinger and C ady 76 and those of Stevens et al 75 As expected by the disagreement in the EOS in Figure 5 our values show disagreement with the bulk modulus reported by both Byrd and Rice 59 and Pastine and Bernecker 77 in Table 3 3.3 UNIAXIAL COMPRESSIONS The princi pal stresses, i.e. the eigenvalues of the stress tensor, as a function of compression ratio are shown in Figure 7 They were arranged such that 1 is the maximum eigenvalue and 3 is the minimum eigenvalue. The large differences in principal stress values between compression directions clearly indicate the anisotropic character of TATB. The compressions along crystallographic directions normal to {100} {010} and {110} do not involve deformations along the c lattice vector, and these directions exhibit greater stresses than compressions which reduce the interplanar distances along the c direction. It is also interesting that the magnitude of 1 i s as high as 43 GPa at a compression ratio of 0.7, twice larger than values calculated for PETN 10 and HMX 73,79 which were approximately 20 GPa.


27 Figure 7 : Principal stresses as a function of V/V 0 The pressure from the hydro static compression calculations is shown for comparison. The animation of the unit cell during compression showed that for the direction normal to {100}, the molecules lose the planar structure at At the gr oups on opposite sides of the molecules bend out of the a b plane. Compressing along the direction normal to {110} we observed the same behavior. The molecules start to deviate from the planar structure at and the NO 2 groups ben d out of the plane at In the direction normal to {010} these two behaviors occur at and respect ively Figure 8 depicts the TATB unit cell at the maximum compression for the direction normal to {110}; one can see deviations from the planar structure and the bending of the NO 2 groups out of the plane. Figure 8 : TATB molecule compressed up to for {110} orientation.


28 Another indication of anisotropic behavior in TATB is observed in the change in energy per atom as a function of v olume. Figure 9 clearly shows that compressions normal to the {100} {010} and {110} planes undergo the greatest change in energy per atom upon uniaxial compression. Qualitatively, this can be explained by strong intramolecular bonds within the a b pla ne as compared to much weaker intermolecular interactions along the c direction. The band gap was also calculated as a function of compression for each direction. As shown in Figure 10 the band gap decreases very slowly as a function of the compression r atio. We compared our band gap calculations with those performed by Wu et al 53 in the direction normal to {001} Wu et al obtained the electronic structure of TATB using the local density approximation (LDA) with norm conserving pseudo potentials. Desp ite the difference in calculation parameters, the results are similar. We found a band gap of 2.35 eV as compared to their band gap of 2.4 eV at zero compression and a band gap of 1.24 eV as compared to 1.4 eV at Figure 9 : Band gap as a function of compression. Figure 10 : Change in energy per atom as a function of compression.


29 The shear stresse s are the driving forces for plastic deformations and mechano c hemical events behind the shock wave front. Therefore, we calculated the maximum shear stresses upon uniaxial compression as where i and j take the values 1, 2, and 3. Th e maximum shear stress max is found by using and but the values of i and j used to find m id and min depend on the principal stresses for the given compression direction. In Figure 11 the maximum shear stre sses as a function of compression are displayed. From the onset of compression up to the shear stresses along the directions normal to {100} {010} and {110} are greater than those for the directions involving compressions along th e c lattice vector. Note the difference in scale for the plots of the shear stresses; max has nearly the same magnitude of m id and more than twice the magnitude of m in Also, TATB has higher shear stresses than both PETN and HMX, with values as high as 18 GPa versus 16 GPa 73 Figure 11 : Shear stresses as a function of compression with max on the left and min on the right. 3.4 CONCLUSIONS First principles simulations of hydrostatic and uniaxial compressions along the crystallographic directions normal to the {001} {010} {011} {100} {101} {110} and {111} planes in TATB were performed. The relaxed structure at zero pressure showed a 1.49% difference in volume as compared with the experimental structure at ambient


30 conditions. The hydrostatic compression data, including pressure, energy per atom, and band gap, were calculated at 0 K and compared with experiment. The calculated bulk modulus is in good agreement with Stevens et al 10% less than that of Olinger and Cady, but shows a large disagreement with the theoretical prediction by Byrd and Rice. Also, our calculations indicate that the a/b ratio remains approximately constant under hydrostatic compression. The uniaxial compression data of this study clearly indicate anisotropic properties in TATB. In particular, the shear stresses for compressi ons normal to {100} {010} and {110} up to are greater than those directions that involve compressions along the c lattice vector. Unfortunately, there are no experimental data regarding the anisotr opic shock sensitivity of TATB. Ho w ever, a recent study 10 by Conroy et al. suggests that there may be a correlation between shear stress and shock sensitivity in EMs. U s ing this assumption we predict that the TATB crystal is more sensitive when compressed normal to the {100} {010} and { 110} planes


31 4.0 HYDROSTATIC AND UNIAXIAL COMPRESSION S OF NITRATE ESTER Re c ently synthesized by Chavez et al ., 80 a novel ener g eti c nitrate ester shows promising properties for use in explosives appli c ations. This material, whi c h we shall refer to as Nitrate Este r 1 (NEST 1) exhibits sensi tivity properties similar to p entaerythritol tetranitrate (PETN), and predi c tions of its performan c e properties indi c ate t hat NEST 1 is as powerful as c y c lotetramethylene tetranitramine (HMX). 80 Further characterizat ion is required to investigate the potential use of this mate rial. Modeling and simulation to understand the properties of NEST 1 will require experimental data that h ave yet to appear in the l iterature. Unfortunately, exper iments designed to measure the isothermal equation of state (EOS), elasti c properties, and sho c k response of an EM can be c ostly to perform. Theoreti c al te c hniques may be applied to predi c t these properties for e ff e c tive c hara c terization of NEST 1 prior to the availability of experimen tal data that c ould take months or years to appear in the literature. Density fun c tional theory (DFT) with van der Waals (vdW) c orre c tions 65 68 c an be used to make reliable pre di c tions for mole c ular c rystals su ch as EMs. In a previ ous w ork by Conroy et al. 60 it was shown that DFT with the empir i c al vdW c orre c tion of Neumann and P errin 67 c an pre di c t the equilibrium volume of several energeti c materials within about 3% error in c omparison with experiment. It was also shown that vdW c orre c ted DFT p rovides bet ter agreement with the experimental EOS at low pressures


32 than DFT alone. We refer the reader to referen c e [67 ] for a detailed des c ription of the empiri c al co rre c tion. In this work, we use vdW c orre cted DFT to investi gate the response of NEST 1 to h ydrostat i c and uniaxial c ompressions perpendi c ular to the {100}, {010}, {001}, {110}, {101}, {011}, and {111} planes. There are two main goals for this study. We w ant to c al c ulate a cc u rate physi c al data for use in larger s c ale simulation of NEST 1 and to invest igate the response of this EM to simulated uniaxial c ompression to investigate the degree of anisotropy in its c onstitutive relationships. We are parti c ularly interested in the shear stresses be c ause similar c al c ulations 10 of the EM PETN indi c ated that th ere might be a c orrelation between shear stresses under uniaxial c ompression and sensitivity to sho c k indu ced deto nation. 4.1 COMPUTATIONAL DETAILS The Vienna Ab initio Simulation Pa c kage (VASP) 43,44 was used to perform all DFT 81 ,8 2 c al c ulations. The PW9 1 fun c tional 69,70 with proje c tor augmented wave (PAW) potentials 71,72 was c hosen be c ause the parameters of the vdW c orre c tion 67 were fi t with these settings. The energy c uto ff for all c al c ulations was 700 eV For ea c h type of c ompression, a Monkhorst Pa c k k grid 74 was used that c orresponded to a spa c ing of 0.08 1 when the unit c ell was c ompressed to half of its experimental v olume. The ele c troni c energy toleran c e was set to 10 6 eV. The c onjugate gradient algorithm was employed to relax all stru c tures until the maximum for c e on any atom in the unit c ell was less than 0.03 eV/. The experimental structure at ambient conditions was relaxed without constraints to determine the calculated equilibrium structure. Hydrostatic compression simulations


33 were per formed by scaling the volume of the unit cell from 100% to 60% of the calculated equilibrium volume. Relaxations of the unit cell shape and atomic coordinates were performed at intermediate steps of 2%. Uniaxial compressions of NEST 1 were performed along seven d irections corresponding the normals of the {100}, {010}, {001}, {110}, {101}, {011}, and {111} planes. For each orientation, the crystal structure was rotated such that the direction of compression was oriented along the x axis. The x component o f each lattice vector was the n scaled by 2% increments from down to At each compression step, the unit cell was held fixed allowing only the atoms to relax. 4.2 EQUILIBRIUM PROPERTIES AND HYDROSTATIC EOS The calculated unit cell properties at equilibrium are displayed and compared wit h experimental values in Table 4 The agreement with the experimental structure is excellent, yielding an error in unit cell volume of approximately 0.1%. Table 4 : Equilibrium unit cell properties of NEST 1. Study a () b () c () Volume ( 3 ) Chavez et al. 8.1228 23.0560 8.5072 90 113.953 90 1456.01 This work 8.119 23.094 8.505 90.00 113.939 90.00 1457.57


34 From the hydrostatic compression simulation, the isothermal equation of state (EOS) was calculated and is shown in Figu re 1 2 with the Birch Murnaghan EOS fit used to obtain the bulk modulus The data fit the EOS very well with the exception of the points below where a change in symmetry of the crystal from to was detected from our calculations at Further work is necessary to determine whether this change in symmetry is an artifact of the vdW correction at high pressure. The lattice constants a b and c as a function of press ure are shown in Figure 13 The lattice constants up to 27 GPa were fit to fifth order polynomials, and the c oefficients are shown in Table 5 Above 30 GPa, each lattice constant shows a change in the trend of the curve, and this corresponds to the change in symmetry predicted at Figure 13 : Lattice parameters of nitrate ester 1 as a function of pressure for the hydrostatic compression simulation. Figure 12 : Isothermal equation of state for nitrate ester 1 (NEST 1). The range in volume is from V/V 0 = 0.60 1.00. The Birch Murnaghan EOS fit used to obtain the bulk modulus is also shown for comparison.


35 Table 5 : Fifth degree polynomial fit s for pressure dependence of lattice constants and unit cell volume. Each constant C n has units of /(GPa) n Parameter C 0 C 1 (x 10 1 ) C 2 (x 10 2 ) C 3 (x 10 3 ) C 4 (x 10 5 ) C 5 (x 10 7 ) A 8.119 1.356 1.455 1.015 3.493 4.551 B 23.09 3.312 3.322 2.174 7.202 9. 149 C 8.505 2.704 3.289 2.285 7.637 9.649 The hydrostatic compression data up to 5 GPa were fit to the Birch Murnaghan 83 Murnaghan, 84 and Vinet 85 equations of state to find the bulk modulus and its pressure derivative. The results for each fit are sho wn in Table 6 Table 6 : Bulk modulus and its pressure derivative for NEST 1 Fitting EOS (GPa) Birch Murnaghan 13.9 8.5 Murnaghan 14.6 6.8 Vinet 14.1 8.1 4.3 UNIAXIAL COMPRESSIONS To investigate the possibility of anisotropic constitutive relationships in NEST 1, uniaxial compressions were performed. For each orientation and uniaxial compression step, the eigenvalues of the stress tensor i.e. the prin cipal stresses, were determine d The greatest of the diagonal elements was taken to be 1 and the lowest 3 The greatest principal stress is approximately the stress tensor element xx which is oriented along the axis of compression. A plot of each of the principal stresses as a fu nction of compression ratio can be seen in Figure 14


36 Figure 14 : Principal stresses 1 (a), 2 (b), and 3 (c) of NEST 1 under uniaxial compression. The principle stresses for the hydrostatic compressi on are included for comparison. The greatest principal stress 1 shown in Figure 14 (a) suggest a low degree of anisotropic sensitivity given the overall, qualitative similarities for each compression orientation. Between and a separation of about 1 2 GPa can be seen between the {011} and {010} orientations, with all ot her orientations lying between. At a deviation from monotonic behavior can be seen for the {010} orientation. Simila rly, other orientations exhibit such deviations at higher compressions. These non monotonic deviations seem to be the result of conformational changes within the unit cell at the compression step of the deviation as seen in animations of the compressions; though it is likely that the atomic configurations that give rise to these deviations are a subset of many possible configurations on the minimized potential energy surface that can result from the conjugate gradient relaxation method. As for the other p rincipal stresses, the change in energy per atom as a function of compression ratio was also calculated for Figure 15 : Change in energ y per atom upon uniaxial and hydrostatic compressions.


37 each orientation as seen in Figure 15 Again, little difference is observed between the orientations studied. The closely distributed curves in bo th Figure 14 (a c) and Figure 15 suggest a relatively isotropic behavior for NEST 1. A comparison of NEST 1 and PETN was made due to the similarities in molecular structure. In contrast to the principal stresses for NEST 1, a greater degree of anisotropy was indicated for the values of the principal stresses with increasing compression ratio in a previous work on PETN under uniaxial compression. 10 Though the two molecules posses similar molecular structure, their crystalline structures are substantially di fferent. Therefore it is not surprising to find that despite the anisotropic behavior of PETN under uniaxial compression, we find NEST 1 to exhibit relatively isotropic behavior compared to PETN. In fact, in comparison with other EMs such as HMX and RDX 79 ,86 we find the same trend. For HMX, the differences are clear when one compares the principal stresses with those found for NEST 1. The spread of the curves for principal stresses and change in energy per atom are notic eably wider for HMX 79 suggesting th at NEST 1 possesses much less of a degree of anisotropic sensitivities in comparison. In the case of RDX however, the principal stresses may not seem to be spread significantly wider than those for NEST 1, but the change in energy per atom does suggest tha t NEST 1 has a lesser degree of anisotropic sensitivity than RDX. Therefore, we predict NEST 1 to be relatively isotropic in shock sensitivity compared to PETN, HMX, and RDX.


38 4.4 CONCLUSIONS First principles simulations of both hydrostatic and uniaxial co mpressions normal to the {100}, {010}, {001}, {110}, {101}, {011}, and {111} planes of the NEST 1 molecular crystal were performed using DFT with the PBE functional, PAW potentials, and an empirical van der Waals correction. The calculated equilibrium str ucture was found to be in excellent agreement with experiment, having an error of about 0.1% for the unit cell volume. The isothermal EOS was calculated and was shown to follow the Birch Murnaghan EOS fit down to At a change in symmetry from to was detected. Using the Birch Murnaghan, Murnaghan, and Vinet EOS fitting schemes for hydrostatic compressions up to 5 GPa, the bulk modulus and its pressure derivative were calculated. The principal stresses and change in energy per atom as functions of compression ratio were calculated for each uniaxial compression in each orientation. Apart from relatively small deviations from monotonic behavio r in the plots of the maximum principal stresses, the principal stresses and energy changes per atom for each of the orientations were more closely distributed than for compression simulations of PETN, HMX, and RDX This suggests relatively isotropic sensi tivities to shock as a fun c tion of compression direction, making NEST 1 unique among other energetic molecular crystals.


39 5.0 HYPERVELOCITY COLLISIONS IN DETONATING PETN An understanding of the mechanisms that bring about the onset of detonation in ene rgetic materials (EMs) has been sought for a number of decades. Originally, chemical decomposition was believed to be more or less describable by classical transition state theory; however, the study of EMs brought a problem to the classical view. When a material undergoes the process of deflagration, the reaction rate of the heated material typica lly increases with temperature. This was found not to be the case for detonatin g EMs; the temperatures associated with a steady state detonation process implied reaction rates that wer e much faster than was observed 12 This finding initiated an exhaustive effort to elucidate the mechanisms responsible for this d e viation from classical behavior seen in detonating EMs. Since the 1970's, efforts have been underway to describe the chemical behavior of nonequilibrium processes such as hyper velocity mo lecular collisions 87 and detonation. 12 17 Decades of work have yielded such elegant concepts as phonon mode c o upling and up pumping 18 21 which provide a basis by which r eaction rates can be predicted for ma n y grained polycrystalline EMs. Unfortunately, these models fall short of explaining the sensitivity anisotropies found in single crystals of EMs, particularly those discovered in PETN by Jerry Dick. 5 The discovery of a nisotropic shock response in EMs brought to light new directions in the study of EM initiation chemistry. Physical models involving hot spot formation began to evolve leading the search for the mechanisms of initiation


40 towards concepts like steric hindran ce along slip systems 6,8,9 or molecular collisions in c r ystal voids and cracks. 2,4,22 Only now, with the rise of computational resources, can such physical models be adequately studied. We present here one such computational study of the hypervelocity co llisions of PETN molecules in order to explore the reaction dynamics that cause the first chemical reactions leading to detonation. In this study, we have determined threshold velocities of initiation for several orientations corresponding to the crystal environment. We have also performed a detailed, qualitative analysis of the role of steric factors and deformations in energy transfer, and calculated the bond dissociation energies (BDEs) required to generate the products observed in simulation. 5.1 COMP UTATIONAL DETAILS The simulation of chemical reactions in the condensed matter phase is a challenging problem because it requires accurate quantum mechanical descriptions of dynamics in a system consisting of a large number of atoms. Therefore, we devise d a model system that contains a manageable number of atoms that are treated fully quantum mechanically using first principles density functional theory (DFT), while preserving essential features of the initiation chemistry. Because the major focus is on chemical dynamics initiated by the shock wave we assume that the initial steps involve the rapid approach and collisions of molecular pairs. Therefore, the pair is extracted from the EM crystal by preserving the local crystalline environment (relative pos ition and orientation). Then, the bimolecular collisions are simulated by assigning the pair of selected molecules the relative velocity and following their collision dynamics.


41 Reactive MD simulations of bimolecular collisions were performed within the Born Oppenheimer approximation using density functional theory (DFT). The trajectory of each atom was calculated by integrating the classical Newton equations of motion: where potential energy and its gradient s are evaluated on the fly by solving the Kohn Sham DFT equations at each time step at fixed nuclear coordinates R i to obtain the total electronic energy Then, which also includes the Coulomb nuclear nuclear repulsion contribution, serves as a potential energy for nuclear dynamics: Equation 1 together with the initial conditions (discussed below) completely specifies the dynamics of the bimolecular collisions. First principles DFT calculations were performed using the linear combination of ato mic orbitals (LCAO) code SIESTA 88 The core electrons are taken into account by employing norm conserving pseudo potentials. The Perdew Burke Ernzerhof (PBE) generalized g radient approximation (GGA) density functional is used. 45 The valence electronic states are expanded using the double zeta plus polarization (DZP) basis set that we specifically optimized in order to achieve the highest level of accuracy in the description of the energies and forces as compared to the fully converged plane wave calculations for several configurations of the colliding PETN molecules. The real space energy cutoff was 200 Ry Equation 3


42 The initial conditions for the reactive collision dynamics simula tions were set up by constructing the two molecule collision complex and specifying initial velocities for each atom of the complex. We are specifically interested in relating the bimolecular reactivity with observed orientation dependent shock sensitivity Therefore, the collisions along the various crystallographic directions within the PETN unit cell were chosen as a starting point for our study. All possible collisions in the following directions were considered: {001}, {010}, {011}, {100}, {101}, {110} and {111}. To simplify this process, each molecule in the crystal lattice was substituted with a sphere of radius R which is half the diameter of the molecule. Figure 1 6 shows this simplified model of the PETN crystal with each molecule indexed for the purpose of identifying a specific collision pair. The inner spheres have a radius of R/ 2 and were devised to distinguish between two different types of collisions: head on and glancing. If the inner spheres of two molecules collide upon translation in a particular direction, then the collision is classified as head on; if only the outer spheres collide, then the collision is labeled as glancing. Once all possible collisions were found, the associated molecules were isolated in their crystallographic ori entation with respect to one another and centered within an empty box. The molecules were then thermally equilibrated to ~300 K Next, one half the Figure 16 : Simplified depiction of PETN lattice using spheres in place of molecules with molecular indices used to designate each collision case. Blue represe nts zones of head on collisions, while grey represents zones of glancing collisions.


43 desired collision velocity was applied to each molecule such that they collided in the center of the box as depicted in Figure 17 Collision velocities ranging from 3 7 km/s were sampled for thirteen different orientations. Due to the computational expense of DFT calculations, a statistical study was not performed. However, the statistics are indirectly samp led by running simulations with multiple collision orientations. 5.2 THRESHOLD VELOCITIES TO INITIATE CHEMICAL REACT IONS The trajectory data obtained from each simulation were analyzed to determine the onset of reaction, the reaction time scale, and the t ype of products. The lowest velocity at which reactions occurred for any particular direction and orientation was taken to be the threshold velocity of initiation for that particular case. All the reactive cases at the threshold velocities share one impor tant feature: the reaction products originate from one molecule while its colliding counterpart r emains chemically intact At higher collision velocities, products were observed to be originating from both molecules. The investigation of the threshold ve locities as a function of initial crystalline orientation allows us to make an important conclusion: the reactive initiation dynamics is orientation dependent. Indeed, we find that each collision orientation has a unique threshold velocity of chemical init iation, see Figure 1 8 Moreover, in addition to having Figure 17 : Collision geometry of two PETN molecules set to collide along a specific direction.


44 different threshold velocities, the collisions along various directions proceed along different reaction pathways as will be made apparent. Collisions normal to the {001}, {010}, and {011} planes ca n be considered to be most sensitive because their threshold velocities are the lowest. The {111} collision case produced reactions at the intermediate threshold velocity of 4.2 km/s Finally, collisions normal to the {010}, {100}, and {101} planes are mo st insensitive with threshold velocities well above 4.0 km/s which is considered to be a typical particle velocity behind the shock wave front. 6 In the study performed by Wu et al 9 comparable threshold velocities were found for the {001}, {100}, and {110} head on collision cases, though glancing collisions showed much lower threshold velocities than those found in our study. This is not unexpected because the glancing collisions performed by Wu et al involved molecules that were initially deformed in orde r to impose external conditions corresponding to slip systems along the maximum resolved sheer stresses As is seen in Figure 1 8 there is no special preference in chemical reactivity for either glancing or head on collisions under ambient conditions. We found a correlation between experimentally observed shock sensitivities reported by Dick, Yoo, et al ., 6,7 and the threshold velocities of reaction initiation, see Figure 18 : Thres hold velocities for reaction initiation as a function of the crystallographic orientation of each collision case. Numbers to the right of each point correspond to the molecular index of specific collisions, see Figure 1 6.


45 Table 7 Essentially, the detonation pressure s shown in the right column of Table 7 are a m easure of the sensitivity for a given direction in bulk PETN. Known sensitive cases with lower detonation pressures have lower threshold velocities than the known insensitive cases with higher detonation pressures. Table 7 : Thresh old velocities, products, times of product formation, and experimentally determined shock sensitivities for bulk PETN. Molecular indices match those in Figure 1 6 Direction Type Molecules Velocity (km/s) N Prod. Products Time to Rx (fs) Exp t. Press (GPa) { 001 } Head On 1 3 3.30 1 NO 2 450 12.1 {010} Glancing 1 2 4.70 1 NO 2 135 N/A Head On 1 4 4.40 1 NO 2 410 { 011 } Glancing 1 2 5.30 1 NO 2 110 N/A 1 3 3.40 1 NO 2 130 1 5 3.20 1 NO 2 225 {100} Head On 1 6 4.70 3 NO 2 ; H 2 CO; NO 2 130; 16 0; 42 0 22.8 { 101 } Glancing 1 2 6.60 4 NO 2 ; NO 2 ; H 2 CO NO 2 90; 100; 316; 480 No Rx P > 31 GPa 1 3 4.60 3 NO 2 ; NO 2 ; H 2 CO 90; 135; 320 1 7 4.40 2 NO 2 ; NO 2 155; 415 { 110 } Head On 1 8 3.90 2 HONO; CO 310; 310 4.2 Glancing 1 2 4.50 1 NO 2 160 { 111 } Head On 1 2 4.20 2 NO 2; NO 2 140; 295 N/A For example, two orientations {001} and {110}, were observed to b e the most sensitive with detonation pressures 12.1 GPa and 4.2 GPa respectively. T he threshold velocities 3.30 km/s and 3.9 km/s for h ead on {001} and {110} collisions are indeed the lowest among all collision cases studied. In contrast, the {100} and {101} case s proved


46 to have high threshold velocity of 4.7 km/s and 4.4 km/s respectively, which correlates with the high detonation pressu re of 22.8 GPa determined by Yoo et al 7 for {100}, and the lack of initiation for {101}. It is important to note that the correlation between detonation pressure and threshold velocity is merely qualitative, especially when comparing the pressures and vel ocities for the sensitive cases which include the so called anomalous elastic respon se of the {110} case at 4 GPa 6 5. 3 PRODUCT FORMATION The number of products formed in reactive threshold cases shown in Table 7 can also serve as an indicator of the di rectional de pendence of chemical initiation as see n in Figure 19 While one might expect more products to form with higher collision velocities, we observed some deviations from this trend. This is evident when comparing {010} mol 1 2 {010} mol 1 4 and { 011} mol 1 2 having threshold velocities of 4.7 4.4 and 5.3 km/s respectively, with {100} {101} mol 1 2 and {101} mol 1 3 having threshold velocities of 4.7 6.4 and 4.6 km/s All co llisions in the first group pro duce only one fragment while all Figure 19 : Number of fragments produced vs. collision veloc ity for hea d on collisions (top) and glancing collisions (bottom ).


47 co l lisions of the second group produce three fragments. In short, the number of products produced by a collision depends as much on the orientation as it does on collision energy. Probing the issue of product formation deeper, one can look at the number of products formed for all collision velocities of any particular direction and see once again that product number does not necessarily scale wit h collision energy, see Figure 19 Though the number of products increases with collision velocity in cases like {111} head on and {110} glancing, other cases such as {100} head on and {101} mol 1 7 glancing show that the number of products can actually diminish at higher collision velocities. One interesting example of this is the {001} head on case in which a singl e product is produced at 3.3 km/s and no products are produced at 3.4 km/s The reason for this is hidden in the details of the dynamics In the 3.3 km/s case, which is discussed later in Figure 21 the reacting nitro group is imparted with translationa l kinetic energy as the co nnecting arm is thrown outward. As the associated N O bond stretches in its oscillation, the connecting formaldehyde group snaps back towards its original configuration, leaving the NO 2 behind. In the 3.4 km/s case, the formaldehy de does not snap back quite as fast due to a twisting deformation in the center of the molecule. As a result, the N O bond that was severed in the 3.3 km/s case is not severed in the 3.4 km/s case. This may be due to the actual positions and velocities o f the individual atoms due to their thermal vibrations. However, the amplitudes of the N O C C and C O thermal vibrations are small (~0.1 ) making it unclear wheth er or not thermal effects are in play. To determine this, more sampling of different equ ilibrated structures is needed.


48 5. 4 REACTION PATHWAYS Several studies have addressed a correlation between bond dissociation energies (BDEs) and impact sensitivities of various collidin g compounds 9 ,89 91 Wu et al conducted such a study for PETN by ca lculating BDEs for PETN using DFT with the BPW91 GGA functional. 9 They found that the N O single bond was weakest with a BDE of 157.8 kJ/mol 9 Another study by Fried et al also found the single N O bond to be weakest with a BDE of 167 kJ/mol 89 We perfor med our own study of the BDEs for PETN including all bonds within the PETN molecular complex. Figure 20 shows these BDEs for the intact molecule, as well as the BDE of the C C bond after the formation of NO 2 As seen in the study of Wu et al ., we found N O 2 production to be the dominant chemical event in our simulations. Every reactive collision performed i n this study began with the di ssociation of the N O bond leading to an NO 2 product as seen in Table 7 T his commonality can be partly attributed to th e relative weakness of the N O bond (182.69 kJ/mol ) in the non reacted PETN molecular complex as seen in Figure 20 It is clear from should be pointed out that some stu dies have suggested that the N O bond is not the first to break when PETN is exposed to a high energy laser, high energy fracture, or high pressure; rather they suggest that the C C and C O bonds are the first to break 93,94 In Figure 20 : Bond dissociation energies for PETN in kJ/mol calculated using PBE functional


49 reference to Ref. [93] whic h cites the high electron density value of the N O bond, we assert that it is the bond order, not the electron density, which deter mines the breaking of the bond. Also, in reference to Ref. [94 ], which finds no peak for NO 2 in a mass spectroscopy study of laser decomposed PETN, we point out that the time resolution required to see the initial stages of chemical initiation is ~1 ps not ~100 The next weakest bond in the PETN complex is the C C bond with a BDE of 250.48 kJ/mol This particular bond scis sion occurred less frequently in threshold cases, but was nonetheless prevalent, occurring only after or during the dissociation of the N O bond in the reacting arm. We calculated the BDE for the C C bond after the dissociation of the accompanying N O bond and found a drop in energy from 250.48 to only 77.07 kJ/mol The lowering of the BDE for the C C bond may explain why H 2 CO forms in the stead of more NO 2 fragments such a reaction pathway is actually favored by way of the BDEs. From Table 7 it is evide nt that some threshold cases produced multiple NO 2 fragments without producing H 2 CO In both such cases for {101} glancing collisions only the outer nitro groups of the colliding PETN molecules interacted leaving the interiors relatively intact. Likewise, the {111} head on threshold collision was oriented such that the nitro groups in multiple arms took the brunt of the collision simultaneously, allowing no close approach to the interior of either molecule. In other collisions in which multiple fragments were formed, close approaches affecting the interiors of the reacting molecule were prevalent. It therefore must be concluded that while BDEs can account for part of the chemical phenomena seen in our simulations, the actual dynamics of the collisions mus t also be considered for a complete picture of the initiation chemistry associated with the detonation of EMs.


50 5. 5 REACTION DYNAMICS study of the dynamic factors of the collisions are necessary to understand the localized transfer of energy leading to the rapid reactions seen in our simulations. As a starting point, we looked at the simple dissociation of the N O bond leading to the formation of NO 2 in single product threshold case s such as the {001} head on collision. The depiction in Figure 21 shows various important events leading to the formation of NO 2 for the {001} collision. Upon collision, the approaching arms in the left molecule are thrown outward Figure 21 (b and c) The nitro group at the top of the left molecule continues to move to the left as the connecting formaldehyde group recoils back to the right Figure 21 (d and e) Because of this opposing motion, the N O bond connecting the nitro group to the formaldehy de group is stretched to the breaking point Figure 21 (f) Essentially, the entire arm was moved out of its equilibrium position, and like a loaded spring tended to return to its original state; however, while the inner part (the formaldehyde group) rec oiled back towards its equilibrium position, the nitro group was prevented from doing so by its momentum. Figure 21 (g) shows a plot of the bond lengths for the N O C C and C O bonds within the reacting arm. The C C and C O bonds remain out of phase duri ng the stretching periods of the N O bond the carbon atom within the arm is moving between the connecting oxygen and center carbon. When the C C and C O bonds finally come into phase, the oxygen atom is pushed to the right as both bonds stretch, severing the N O bond. Then the carbon atom is pulled back down as both the C C and C O bonds contract, making the N O scission permanent.


51 Figure 21 : {001} head on mol 1 3, 3.3 km/s : (a f) Collision dynamics leading to an N O bond sci ssion; (g) bond lengths plotted against time; vertical lines show frames depicted in (a f).


52 In general, the loading of the bonds in the reacting molecule is mainly localized to only those parts immediately affected by the collision. A prime example of th is is the {100} head on threshold case, shown in Figure 22 in which an N O bond is immediately broken as a nitro group from the right molecule strikes the reacting nitro group in the left molecule almost laterally. T he N O bond scission is immediately fol lowed by a C C bond scission and subsequently by another N O bond. The N O bond scission happens fast, occurring before any massive deformation of the rest of the reacting molecule Figure 22 (a c) The stretching mode of the N O bond in the non reacti ng molecule on the right is loaded as it pushes against the first reacting nitro group on the left. As seen in Figure 22 (g) the amplitude of the N O stretching mode in the non reacting molecule increases once the N O bond in the reacting molecule is clea ved. With the majority of the collision energy localized in the reacting arm, the connecting formaldehyde group is pulled down, stretching the associated C C bond to its breaking point Figure 22 (d) As discussed earlier, The BDE for the N O bond is ini tially the lowest in the molecular complex; however, once the N O bond breaks, the BDE for the associated C C bond drops significantly to become the lowest (from 250.48 down to 77.07 kJ/mol ) which can in part account for the reaction pathway seen here. I n Figure 22 (e) the upper arm of the reacting molecule begins a downward swing, bending sharply between the central Carbon atom and the formaldehyde group. The bend is so sharp that the Oxygen atom in the formaldehyde group bonds with the central Carbon at om while the downward motion of the nitro group and the stretching of the associated N O bond leads to the formation of another NO 2 Figure 22 (f) The fact that this reaction pathway occurred as described offers supporting evidence that the chemical


53 Fi gure 22 : {100} head on mol 1 6, 4.7 km/s : (a f) Collision dynamics leading to N O, C C, and N O bond scissions; (g) bond lengths plotted against time; bond length of suppressed, non interacting N O bond responsible for first react ion is shown by dotted line.


54 initiation due to a hypervelocity collision is not thermal. If the reactions were thermal, then the formation of NO 2 and H 2 CO would not have been immediate; at the very least, they would have had to occur only after some numbe r of oscillations of their associated, excited modes of vibration. It is also important to note that between frames (d) and (f) in Figure 22 some of the nitro groups in the non reacting molecule seemed to dissociate temporarily as seen in Figure 22 (f) If these apparent dissociations, which were identified by extremely large bond lengths, can indeed be considered reactions, then these nitro groups reform as they travel along trajectories similar to their parent molecule. The rarest reaction pathway ob served in the threshold collision cases involves C H bond scissions at a rather low collision velocity. In the {110} head on threshold case, NO 2 is formed in much the same manner as in the {001} threshold case. H 2 CO forms almost immediately after the form ation of NO 2 in a manner similar to that of the {100} threshold case. These products form almost simultaneously and have rotations imparted to them. As a result, the h ydrogen atoms in the H 2 CO come within bonding distance of both the central Carbon atom o f the reacting molecule, and an Oxygen atom in the NO 2 causing the bond scissions as depicted in Figure 23 (a d). The C H reaction events can be seen in Figure 23 (e) where the plots of the bond lengths of the reactants cross those of the resulting prod ucts. The thermal behavior of the hydrogen bonds is evident before and after the reaction events. It is interesting to point out that a study done by Dick, which reported an anomalous shock response for the {110} direction at 4.2 GPa 6 also reported Raman spectra data that suggested the formation and decomposition of HONO The {110} head on collisions were the only set


55 of collision that produced HONO with collision velocities below 5 km/s with the {111} head on collisions being the only other set to produ ce HONO below 6 km/s Because C H thermal vibrations were at play in the formation of CO and HONO and because of the unique reaction dynamics, we performed another simulation of the {110} head on collision at the threshold velocity of 3.9 km/s In this simulation, we equilibrated the bimolecular complex to 310 K in order to investigate the possibility that the reaction Figure 23 : {110} head on mol 1 8, 3.9 km/s : (a d) Collision dynamics leading to C H reaction event; (e) bond lengths of both reactants and products are plotted over entire time span.


56 pathway we observed in the first simulation may have been a statistical aberration. We found that the sequence of chemical events in the 310 K case were nearly identical to those in the 300 K case, differing only in that the events occurred about 15 fs earlier in the later simulation. This suggests that thermal factors leading up to the reaction within a comparable temperature range do not affect the reaction pathway in this case t he formation of HONO and CO is an intrinsic property of the {110} head on collision case at the threshold velocity. The C H scissions occurred due to a combination of the relatively slow rotation of the H 2 CO the fast thermal vibrations of the C H bonds, a nd the proximity of the H 2 CO to the NO 2 and parent molecule. When we looked at similar reactions involving C H scissions in some of the high energy collisions, we found the reaction mechanism to be roughly the same either an oxygen atom from a fragment or a central carbon from a reacted PETN molecule passes within bonding distance of a hydrogen atom and steals it away from the connecting carbon in a formaldehyde group. 5. 6 NON EQUILIBRIUM CHEMISTRY Naturally, if it is assumed that the reactions caused from hypervelocity collisions of PETN are due to a combination of mechanical deformation, steric proximity, and localized vibrations (excited or thermal), then classical transition state theory should not hold. In Table 7 we see that reactions occur betw een 100 fs and 500 fs which puts the reaction timescale ~10 13 s exactly what was sug gested by Walker et al 16,17 This is important since it is known from Arrhenius calculations that thermal equilibration should take 10 11 s If one were to assume that t he reactions are temperature driven, and consider a bimolecular collision of PETN at 4 km/s in which the energy is distributed throughout


57 all of the internal degrees of freedom within the binary collision complex, then the system would have a temperature o f 2000 K The reaction rate of the N O bond scission which produces the NO 2 products could then be calculated as follows: Taking the N O stretching mode frequency as and the bond dissociation energy as kJ/mol the rate of reaction would be or one reaction every ns This is much too slow for the initial chemical events leading to detonation since the time to detonation for PETN is 6 which would only allow for about 154 reactions to occur during the run to detonation. In addition, this reaction timescale is several orders of magnitude larger than those seen in simulation. Finally, the orientation dependence of these co llisions strongly suggests that thermal processes are not at work here. If the reactions were temperature driven, then energy would be evenly distributed throughout the colliding molecules, thus any difference in threshold velocities for all possible colli sion orientations would be statistical in nature. But, the correlation between threshold velocities and detonation pressures suggests that the orientation dependence is real and not a statistical aberration. From this it must be concluded that these initia l chemical events are not temperature driven; rather they are driven directly by the physical dynamics of the collisions due to orientation in a hypervelocity chemical regime. E quation 4


58 5. 7 CONCLUSIONS First principles reactive molecular dynamics simulations of b imolecular PETN collisions revealed that the sensitivities of the various collision orientations given by the threshold velocities of initiation correlate to the experimental sensitivities determined for bulk crystalline PETN. 5,7 Collisions normal to the { 001}, {011}, and {110} planes are most sensitive with lowest threshold velocities of 3.2 3.3 and 3.9 km/s r espectively, while collisions normal to the {010}, {101}, and {100} planes are most insensitive with threshold velocities of 4.4 4.4 and 4.7 km/s respectively. The {111} collision case is at the middle of the spectrum with an intermediate threshold velocity of 4.2 km/s The reactive collisions suggest the formation of NO 2 as the dominant reaction pathway in all cases, with H 2 CO formation for cases in which a large amount of collision energy is localized in a particular arm. In such cases, we find th at the BDE for the C C bond drops after the adjacent NO 2 is liberated. T he sensitive {110} mol 1 8 threshold case produced both HONO and CO due to bot h the proximity of the hydrogen atoms in the H 2 CO product to the NO 2 and parent molecule as well as C H thermal vibrations. A relationship between the steric orientation of the colliding molecules and c hemical enhancement is evident by both threshold velo cities of initiation and product formation. Reactive cases yielded product formation on timescales of the same order of magnitude as the periods of vibration modes with in the PETN molecule, ~10 13 s These reaction timescales are much too short to be drive n by thermal processes, thus leading us to conclude that the direct physical dynamics of the collisions are paramount to initiation.


59 6.0 SUMMARY 6.1 SHOCK COMPRESSION OF DIAMOND DFT s hock compression calculations of diamond in the <100>, <110>, and <1 11> directions were performed both to test the accuracy of the REBO potential under high compression, and to determine the validity of the anomalous shock response found in MD simulations 34 of shocked diamond using REBO. For compressions in <110> and <111> we see values for the shear stresses approach zero and then increase as the compression ratio is increased. A shockwave of appropriate strength to produce such a compression propagating through the diamond lattice would cause an elastic response despite the fact that shocks of lesser strength caused plastic deformations. This is the anomalo us e lastic response which was observed in a previous MD study 34 of shocked diamond using the REBO potential which showed stable elastic behavior in the same compression region as seen where the shear stresses approach zero for the DFT data in Figure 3. Though it is likely that such an anomalous elastic regime would be followed by plastic deformations behind the shockwave front, it is equally likely that the anomalous elastic response is an observable phenomenon given the stability of the MD calculations. 6.2 SHOCK COMPRESSION OF TATB DFT calculations of the equilibrium properties, hydrostatic EOS, and bulk modulus were performed for TATB and found to be in excell ent agreement with


60 experiment. 75 77 This agreement was contingent upon the use of the empirical va n der Waal correction with DFT. 60 Uniaxial compressions were then performed by first rotating the crystal such that the direction of compression was aligned with the x axis, and then rescaling the lattice vectors in the x direction. We found highly anisotr opic behavior for TATB upon uniaxial compression, with higher values for the principal stresses, shear stresses, and energy change per atom for the directions normal to the {100}, {010}, and {110} planes. It is interesting to note that these particular ori entations involve no compressions along the c lattice vector. Through the assumption of positive correlation between stress and sensitivity, we predict that TATB is most sensitive to compressions normal to the {100}, {010}, {110} planes. 6.3 SHOCK COMPRESS ION OF NITRATE ESTER To date, very little experimental or theoretical work has been done for the newly synthesized nitrate ester 1, or NEST 1. 80 For this reason we performed DFT calculations with an empirical van der Waals correction to predict its mechan ical properties. We first performed hydrostatic compressions of the bulk crystal and fit the data to the Birch Murnaghan EOS 83 to determine the bulk modulus. Then we performed uniaxial compressions in the directions normal to the {001}, {010}, {011}, {100 }, { 101}, {110}, and {111} planes. From the uniaxial compression data, we calculated the principal stresses and energy change per atom as a function of compression ratio We found very little difference between the values of each co mpression direction indicating relatively isotropic behavior. Assuming that stress and sensitivity are correlated, we predict that NEST 1 has no enhanced sensitivity in any of the directions we investigated.


61 6.4 HYPERVELOCITY COLLISIONS IN DETONATING PETN We performed first principles MD simulations of colliding PETN molecules over a 0.5 ps time frame to elucidate the mechanisms of chemical initia tion that leads to detonation. From our collision simulations we catalogued the threshold velocities of initiat ion, reaction times, and product formation for each case We found enhanced chemical sensitivity for collisions corresponding to the crystallographic directions normal to the {001}, {011}, and {110} planes, and high insensitivity for directions normal to t he {010}, {100}, and {101} planes. The sensi tivities found for the {001}, {110}, {100} and {101} orientations correlate to known experimental sensitivities. 5 7 We found the production of NO 2 to be prevalent in all reactive cases, with H 2 CO the second most prevalent reaction. A study of the bond dissociation energies of PETN revealed that the N O single bond which breaks to form NO 2 has the lowest BDE in the unreacted PETN molecular complex with a value of 182.69 kJ/mol T he bond dissociation energy of the C C bond connecting the formaldehyde group to the cen ter carbon atom has a value of 250.5 kJ/mol for the unreacted molecule, but reduces to 77.1 kJ/mol after the adjacent NO 2 group on the same arm of the PETN molecule is liberated by N O bond scissi on. Th is finding suggests that once NO 2 is formed, the preferred reaction pathway becomes the formation of H 2 CO. In addition to the formation of NO 2 and H 2 CO, we found a third reaction pathway for the threshold case of {110} which evolves from the initial form ation of the aforementioned pr oducts. This pathway involves both the thermal vibrations of the hydrogen atoms on the formaldehyde group and the proximity of these atoms to the


62 parent PETN molecule and the NO 2 product. This is particularly interesting cons idering that the BDE of the C H bond is very high ~285.4 kJ/mol This phenomenon is not a statistical aberration as it occurred for two differently equilibrated runs, and since HONO has been observed in spectroscopic studies 6 of PETN detonated by shock ap plied in the <110> direction. Finally, our reactive simulations of colliding PETN molecules validate the assumption that thermal decomposition does not occu r in initiation. We find that the reaction timescales are ~10 13 s as predicted by various theoret ical works. 16,17 Because thermal equilibration requires timescales of ~10 11 s we find that it is the direct dynamic al interactions of the colliding molecules that are responsible for initiation.


63 LIST OF REFERENCES 1. R. Courant and K.O. Friedrichs, Sup ersonic Flow and Shock Waves Interscience, New York, 1948. 2. F.P. Bowden, K. Singh, Nature 172, 378, 1953. 3. E.K. Rideal and A.J.B. Robertson, Proc. Roy. Soc., A, 195, 135 (1948). 4. B.L. Holian, T.C. Germann, J.B. Maillet, C.T. White, Phys. Rev. Lett. 89.285501 (2002) 5. J.J. Di c k, Appl. Phys. Lett. 44, 859, (1984) 6. J.J. Di c k, R.N. Mulford, W.J. Spen c er, D.R. Pettit, E. Gar c ia, D.C. Shaw, J. Appl. Phys. 70 (7), 3572, (1991) 7. C.S. Yoo, N.C. Holmes, P.C. Sours, C.J. Wu, F.H. Ree, J. Appl. Phys. 88 (1), 70, (2000) 8. J .J. Di c k, Sho c k Initi c ation Sensitivity of PETN: A Steri c Hindran ce Model, LA UR_91 3413 (1991) 9. C.J. Wu, F.H. Ree, C.S. Yoo, Propellants, Explosives, Pyrote c hni c s 29 (2004), No. 5 10. M.W. Conroy, I.I. Oleynik, S.V. Zybin, and C.T. White, Phys. Rev. B 77, 094 107 (2008). 11. G.A. Peterson, J. Chem. Phys., 3, 107. 12. H. Eyring, S c ien c e, 199 (4330), 740, (1978) 13. H. Eyring, F.E. Walker, S.M. Ma, N. Coon, Pro c Natl. Ac ad. S c i. USA, 77, 2358 (1980) 14. C.S. Coff ey, E.T. Toton, J. Chem. Phys., 76 (2), 949, (1982) 15. F.J. Zerill i, E.T. Toton, Phys. Rev. B, 29 (10), 5891, (1984)


64 16. F. Walker, J. Appl. Phys. 63 (11), 5548, (1988) 17. F.Walker, Propellants, Explosives, Pyrote c hni c s 19, 315 326 (1994) 18. D.D. Dlott, J. Opt. Soc. Am. B Vol. 7, No. 8 (1990) 19. H. Kim, D.D. Dlott, J. Chem. Phys. 93 (3), (1990) 20. D.D. Dlott, J. Opt. So c Am. B 7, No. 8, (1990) 21. D.D. Dlott, Annu. Rev. Phys. Chem. 50, 251 278, (1999) 22. tudes de Dynamique Chimique 1884. 23. Properties of Natural and Synthetic Diamond edited by J.E. Field (Academic Press, 1 992). 24. H.K. Mao and R.J. Hemley, Nature 351, 721 (1991). 25. A.L. Ruoff H. Luo, and Y.K. Vohra, J. Appl. Phys. 69, 6413 (1991). 26. See for example, High Pressure Shock Compression of Solids volumes I III, edited by L. Davison, and M. Shahinpoor (Springer, 1998 ). 27. J.B. Donnet, E. Fousson, T.K. Wang, M. Samirant, C. Baras, M. Pontier Johnson, D iamond and Related Materials 9, 887 (2000). 28. G.R. Willmott, W.G. Proud, and J.E. Field, J. Phys. IV France 110, 833 (2003). 29. D.K Bradley, J.H. Eggert, D.G. Hicks, P.M. Cellie rs, S.J. Moon, R.C. Cauble, and G.W. Collins, Phys. Rev. Lett. 93, 195506 (2004) 30. H. Nagao, K. G. Nakamura, K. Kondo, N. Ozaki, K. Takamatsu, T. Ono, T. Shiota, D. Ichinose, K. A. Tanaka, K. Wakabayashi, K. Okada, and M. Yoshida, Phys. Plasmas 13, 052705 ( 2006) 31. S. Brygoo, E. Henry, P. Loubeyre, J. Eggert, M. Koenig, B. Loupias, A. Benuzzi Mounaix, and M.R. Le Gloahec, Nature Materials, 6, 274 277 (2007). 32. D.W. Brenner, O.A. Shenderova, J.A. Harrison, S.J. Stuart, B. Ni, and S.B. Sinnott, J. Phys. Cond. Ma tt. 14, 783 (2002). 33. D.W. Brenner, Phys. Rev. B 42, 9458 (1990). 34. K. McLaughlin, I.I. Oleynik, S.V. Zybin, M.L. Elert, and C.T. White, AIP Conference Proceedings 955, 321 (2007).


65 35. S.V. Zybin, M.L. Elert, and C.T. White, Phys. Rev. B 66, 220102 (2002). 36. S.V. Zy bin, I.I. Oleynik, M.L. Elert, and C.T. White, in MRS Proceedings "Synthesis, Characterization and Properties of Energetic/Reactive Nanomaterials", 800, AA7.7 (2003). 37. P.E. van Camp, V.E. van Doren, and J.T. Devreese, Solid State Comm. 84, 731 ( 1992). 38. J.J. Zhao, S. Scandolo, J. Kohano_, G.L. Chiarotti, and E. Tosatti, Appl. Phys. Lett. 75, 487 (1999). 39. R.H. Telling, C.J. Pickard, M.C. Payne, and J.E. Field, Phys. Rev. Lett. 84, 5160 (2000). 40. H. Chacham, and L. Kleinman, Phys. Rev. Lett. 85, 4904 (2000). 41. D. Rou ndy, and M.L. Cohen, Phys. Rev. B 64, 212103 (2001). 42. O.H. Nielsen, Phys. Rev. B 34, 5808 (1986). 43. G. Kresse, and J. Furthmuller, Phys. Rev. B 54, 11169 (1996). 44. G. Kres se and J. Furthmuller, Comput. M ater. Sci. 6, 15 (1996). 45. J.P. Perdew, K. Burke, and M. Ern zerhof, Phys. Rev. Lett. 77, 3865 (1996). 46. K. Kunc, I. Loa, and K. Syassen, Phys. Rev. B 68, 094107 (2003). 47. D.C. Swift, and G.J. Ackland, Appl. Phys. Lett. 83, 1151 (2003). 48. B.M. Dobratz, The Insensitive High Explosive Triaminotrinitrobenzene (TATB): Develop ment and Characterization 1888 1994, Los Alamos National Laboratory, LA 13014 H, August 1995. 49. S.F. Rice, R.L. Simpson, The Unusual Stability of TATB: A Review of the Scientific Literature, Lawrence Livermore National Laboratory, Livermore, CA, Report. 50. H. C ady and A.C. Larson, Acta Crystallogr.18, 485 (1965). 51. M. Pravica, B. Yulga, Z. Liu, and O. Tschauner, Physical Review B 76 064102 (2007). 52. M. Riad Manaa, Appl. Phys. Lett. 83, 1352 (2003).


66 53. C.J. Wu, L.H.Yang, and L. E. Fried, Physical Review B, 67, 235101 (2003). 54. J.J. Dick and J.P. Ritchie, J. Appl. Phys. 76, 2726 (1994). 55. J.J. Dick, J. Appl. Phys. 81, 601 (1997). 56. F.Williams, Adv.Chem.Phys. 21 ,289 (1971). 57. J.J. Gilman, Philos. Mag. B 67,207 (1991). 58. E.F.C. Byrd, G.E. Scuseria, and C.F. Chabalowski, J.Phys.Chem. B108,13100 (2004). 59. E.F.C Byrd and B.M. Rice J.Phys.Chem. C111,2787 (2007). 60. M. W. Conroy, Y. Lin, I. I. Oleynik, C. T. White, submitted for publication. 61. J. F. Dobson, in Topics in Condensed Matter Physics, edited by M. P. Das (Nova, New York, 1994). 62. W. Kohn, Y. Meir, and D. E. Makarov, Phys. Rev. Lett. 80, 4153 (1998). 63. K. Rapcewicz and N. W. Ashcroft, Phys. Rev. B 44, 4032 (1991). 64. Y. Andersson, D. C. Langreth, and B. I. Lundqvist, Phys. Rev. Lett. 76, 102 (1996). 65. M. Els tner, P. Hobza, T. Fraunheim, S. Suhai, and E. Kaxiras, J. Chem. Phys. 114, 5149 (2001). 66. S. Grimme, J. Comput. Chem. 25, 1463 (2004). 67. M. A. Neumann and M. A. Perrin, J. Phys. Chem. B 109, 15531 (2005). 68. Q. Wu and W. Yang, J. Chem. Phys. 116, 515 (2002). 69. J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh, and C. Fiolhais, Phys. Rev. B 46, 6671 (1992). 70. J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh, and C. Fiolhais, Phys. Rev. B 48, 4978 ( 1993). 71. P.E. Blochl, Phys. Rev. B 50, 17953 (1994). 72. G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).


67 73. Conroy, M. and Oleynik, I. I. and Zybin, S. V. and White, C. T. (2007) Anisotropic constitutive relationships in energetic materials: PETN and HMX. I n: Shock Compression of Condensed Matter 2007: Proceedings of the Conference of the American Physical Society Topical Group on Shock Compression of Condensed Matter, Waikoloa, HI, 2429 June 2007. American Institute of Physics Conference Proceedings (955 ). American Institute of Physics, Melville, NY, pp. 361 364. 74. H.J. Monkhorst and J.D. Pack, Phys. Rev. B 13, 5188 (1976). 75. L. L. Stevens, N. Velisavljevic, D. E. Hooks, D. M. Dattelbaum, Prop., Explos., Pyrotech. 4/2008,(p 286 295). 76. Olinger, B.W.;Cady, H.H. The Hydrostatic Compression of Explosives and Detonation Products to 10 GPA (100 kbar) and their Calculated Shock Compression: Results for PETN, TATB, CO2 and H2O;Sixth Symposium on Detonation, San Diego, CA,Aug 24 27,1976; Los Alamos LA UR 76 1174. 77. D.J P astine, R.R Bernecker, P,V,E,T equation of state for 1,3,5 triamino 2,4,6 trinitrobenzene. J.Appl. Phys., 45:4458 4468,1974. 78. C.K. Gan, T.D. Sewel, and M. Challacombem, Phys. Rev. B 69, 035116 (2004). 79. M.W. Conroy, I.I. Oleynik, S.V. Zybin, and C.T. White, J Appl. Phys. 104, 053506 (2008). 80. D. E. Chavez, M. A. Hiskey, D. L. Naud, and D. Parrish, Angew. Chem. Int. Ed. 47, 8307 (2008). 81. P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964). 82. W. Kohn and L.J. Sham, Phys. Rev. 140, A1133 (1965). 83. F. Birch, Phys. Re v.71, 809 (1947). 84. F. D. Murnaghan, 'The Compressibility of Media under Extreme Pressures', in Proceedings of the National Academy of Sciences 30, 244 (1944). 85. P. Vinet, J. Ferrante, J.R. Smith, and J.H. Rose, J. Phys. C 19, L467 (1986). 86. M.W. Conroy, I.I. Ol eynik, S.V. Zybin, C.T. White, J. Appl. Phys., 104, 113501 (2008). 87. J. Durbin, Annu. Rev. Phys. Chem. 1973.24:97 120


68 88. Spanish Initiative for Thousands of Atoms, 89. L.E. Fried, M.R. Manaa, P.F. Pagoria, R.L. Simpson, Annu. Rev. Mater. Res. (2001) 31:291 321 90. C. Cao, S. Gao, J. Phys. Chem. B (2007), 111, 12399 12402 91. M. Fang, L. Zhe, Y. Fu, Chinese Journal of Chemistry (2008) 26, 1122 1128 92. E.A. Zhurova, A.I. Stash, V.G. Tsirelson, V.V. Zhurov, E.V. Bartashevich, V.A. Potemkin, A.A. Pinkerton, J. Am. Chem. Soc., 2006, 128 (45), 14728 93. W.L. Ng, J.E. Field, H.M. Hauser, J. Appl. Phys. 1986, 59, 3945 94. Jmol, 95. RasMol,


ABOUT THE AUTHOR Aaron C. Landerville was awarded a B.S. in physics from the University of South Florida (USF) in May of 2007. One year before graduating, he joined the Materials Simulation Laboratory (MSL) at USF with support from the NSF REU fellowship, and began work on computa tional modeling of energetic materials (EM s ) Since joining MSL work has been presented at the annual March meeting of the American Physical Society (APS) in 2007 and 2008, at the bi annual APS Shock Compression of Condensed Matter meeting in June, 2007 and at the Frontiers in Nanoenergetics AFRL workshop, Eglin AFB, October, 2008 Aaron also presented at the USF Undergraduate Research Symposium in 2007 where he received an award for best presentation in the physical sciences and at the Annual USF Graduate Symposium in 2007 and 2008. Aaron is currently pursuing his doctorate at USF under the supervision of Ivan I. Oleynik.