Predicting the reactivity of energetic materials: an ab initio multi-phonon approach

The impact sensitivity of energetic materials is successfully predicted using an ab initio model based on the concepts of phonon up-pumping.


Introduction
Energetic materials (EMs; explosives, propellants and pyrotechnics) are characterized by their ability to convert large amounts of chemical potential energy rapidly into kinetic energy. This has made EMs central to the development of many technologies for both civilian and military application. Paramount to these technologies is the safe and controlled initiation of the EM, with accidental initiation being potentially life-threatening. 1 The need for stringent safety and performance controls has greatly limited the number of compounds used in EMs. While many new EMs have been proposed, 2 successful developments typically favour the preparation of composite materials that consist of a mixture of energetic compounds and additives. Amongst others, this approach has culminated in the formulation of dry mixtures, or more commonly, polymer-bonded explosives (PBXs) such as Semtex and C4. Recent paradigm changes in EM research have motivated focused studies towards development of insensitive munitions that exhibit minimal environmental impact. 3 The development of new energetic molecules has been paramount to these efforts, accompanied by growing interest in the development of multi-component energetic solids. 4 The study of many important physical properties of EMsincluding detonation velocity, enthalpy of formation, heat of explosion, and temperature of detonation -is amenable to computational approaches and quantitative predictions. [5][6][7] The ability to predict properties quantitatively allows for powerful strategies to screen for new energetic molecules. 8 Until recently, and despite a thorough understanding of materials chemistry and detonation physics, there remained very little understanding of the EM initiation mechanism. Therefore, the sensitivity of EMs could not be predicted reliably. Early work by Bowden and Yoffe 9,10 showed that the bulk heating associated with mechanical impact is insufficient to describe the initiation processes. Thus, the concept of 'hot-spots' was developed to describe the localization of energy into highly concentrated regions capable of inducing bond dissociation. A number of more sophisticated mechanisms have since been proposed that typically describe the concentration of energy at interfaces (e.g. grain or particle boundaries) and the generation of high local temperatures or stresses. The former is particularly problematic in the description of processes in molecular materials, where melting occurs at temperatures considerably lower than those required to induce bond rupture. 11 Regardless of the hot-spot origin, the mechanism by which energy leads to bond rupture has yet to be elucidated.
Considerable effort has been devoted to identifying a set of parameters capable of rationalizing the impact sensitivity of EMs. 12 Noting that EM initiation necessarily involves bond rupture, models that are based on isolated molecules have often aimed at correlating impact sensitivity with the stability of bonds in the explosive molecule. This is done, for instance, by computing dissociation energies 13 or analysing the topology of the electron density. 14 However, these models are incapable of describing correctly the different sensitivity of polymorphic forms, i.e. crystalline phases that differ in the spatial arrangement of the energetic molecules. For example, -HMX is much more sensitive to impact than -HMX. 15 Furthermore, the impact sensitivity of an EM can be widely tuned by multi-component crystallization. 4 In such systems, the structure of the energetic molecule itself again remains unaltered. This problem highlights that improved models of EM initiation must also account for characteristic parameters pertinent to the solid state. A number of such parameters have been suggested, including crystal void space and compressibility, 16 intermolecular distances, 17 and intermolecular interaction energies. 18 In models that take account of these characteristics, the process of initial bond rupture has been identified as central. Authors have thus attempted to correlate bond rupture propensity to the electronic band gap 19 and, in turn, impact sensitivity. However, as these methods offer only a qualitative rationale for impact sensitivity, they are generally restricted to chemically similar energetic molecules. In view of this restriction, other attempts have included a consideration of the dynamics of a mechanical impact, for example, by studying the effects of stress on band gaps and bond dissociation energies, [20][21][22] or by considering the interlayer energies associated with imposing a structural distortion within the crystal. 23,24 These more recent methods have typically been restricted to very small sample sets consisting of only a few EMs.
Coffey and Totton, 25 and subsequently Zerilli and Totton, 26 described a mechanism for the localisation of energy in crystalline materials that is based on vibrational energy transfer. This model assumes that mechanical energy, which excites low frequency phonon modes, undergoes rapid up-conversion ('up-pumping') through a process of phonon-phonon collisions. This ultimately increases the population (and hence amplitude) of internal molecular modes beyond covalent bond dissociation limits. The model was then developed further by Dlott and co-workers, 27 who included a description of compression 28 and localization of energy at crystal defects. 29 Various authors have also developed predictive models of impact sensitivity which are based on the phonon uppumping description of mechanically perturbed solids. The earliest attempts of this type used data from inelastic neutron scattering spectra, 30,31 and were followed by work based on experimental Raman spectra. 32 Some work was based on an analysis that considered the vibrational multi-phonon density of states, 30,33,34 while other attempts have instead focused on up-pumping by highorder harmonic (e.g. overtone) coupling processes. 32 Despite the promise offered by a model rooted in multi-phonon up-pumping to rationalize impact sensitivity of EMs, models based on experimental data alone do not offer predictive value. Correspondingly, the first results of vibrational up-pumping methods that rely on an ab initio input based on gas phase 35 and periodic structure models 34,36,37 have now appeared. Using zonecentre vibrational frequencies for molecular crystals, Bernstein 36 was able to rank the relative impact sensitivity of six chemically similar organic EMs. In our previous work, 37 we have investigated the vibrational up-pumping in nine azide-based crystals by calculating the complete vibrational dispersion curves. A correct ranking of these materials was achieved when both overtone and combination pathways were taken into account in the up-pumping process. Thus, while certain vibrational up-pumping models appear to be capable of describing impact sensitivity in EMs, it is not yet clear which models are applicable to a more chemically diverse set of complex organic EMs. This problem motivates the work in the present manuscript that sets out to develop further an ab initio model to describe the impact sensitivity for a set of chemically diverse EMs.

General Theory
In this paper, the term phonon is used in reference to the delocalised, external vibrational modes of a crystal. Vibrational modes that are localised on molecules are denoted as internal vibrational modes.
A mechanical perturbation can lead to the excitation of a crystalline material through various mechanisms. A high-energy impact, for instance, will lead to compression of a material and associated vibrational mode hardening. The magnitude of hardening depends on the average material Grüneisen parameter, Γ = (1/ ) , where and denote frequency and volume, ∑dln(ω)/dln( ) ω respectively, and N is the number of vibrational modes. Within the harmonic approximation, a material with energy = 0 + ℏω will have a heat capacity at constant volume, ( + 1/2) , and is the phonon density of β = ( ) -1 ϵ = ℏω (ω) states. This rationalises the simultaneous increase in and decrease in , which necessitates an increase in the temperature of the material, . A similar logic applies to the Maxwell relation, in which is the isothermal compressibility that relates to the internal energy via the relationship . Here is ∂ derived from an increase in the temperature-dependent phonon pressure. Hence Equation 2 can also be regarded as a description of an increase in temperature that is proportional to . This / description leads to adiabatic heating of the bulk material according to: 38 Equation 3.
It follows that highly compressible materials undergo greater adiabatic heating. The Grüneisen parameters tend to be considerably larger for external vibrational modes than for internal states, and hence the initial adiabatic heating occurs predominantly in the external modes. As the heat capacity of the phonon modes is markedly lower than the total heat capacity, the phonon quasitemperature greatly exceeds that of the final bulk temperature given by Equation 3.
Following initial excitation, the excess phonon energy dissipates through phonon-phonon collisions. By truncating the mathematical description of this process at the leading anharmonic term, the three-phonon collisions (i.e. the combination of two phonon states to create a third) dominate. This can be described by the general Hamiltonian The third term of Equation 4 refers to the conversion of energy between phonon states (ph) and internal vibrational states (v), and is generally defined by a Golden rule formalism of the form, in which the bracketed term represents the two phonon density of states and offers a count for the number of coupling pathways that link a high-frequency mode ( to two lower frequency modes ( Ω) 1 and is the cubic anharmonic coupling coefficient, which 2 ). (3) describes the strength of interaction between three vibrational modes. For the coupling of three external modes, q, It should be noted that Equation 6 can also be formulated for the coupling of internal vibrational modes (Q), or a mixture of external and internal modes. A number of different approaches have been adopted to approximate , 33,34 including its explicit (3) calculation. 39 Unfortunately, the latter is unfeasibly expensive for large systems. A common approach to reduce computational cost is the use of the average anharmonic approximation, which is valid since a large number of coupling states exist within a solid material. 40 Since combinations of modes with anomalously large solutions of Equation 6 are offset by those with anomalously small solutions, can be regarded as representative of the overall 〈 (3) 〉 coupling strength. Furthermore, an explicit consideration of symmetry is unnecessary because the majority of scattering processes do not occur at the Brillouin zone centre or boundaries. Therefore, it has been suggested that the magnitude of can be approximated by taking account of the type of vibrational modes involved. 40 It is reasonable to assume that external modes (q) scatter more efficiently than internal modes (Q) and therefore the magnitude of decreases in the series as qqq > qqQ > qQQ It is generally found that the addition of a fourth phonon (i.e. qqqq) is equivalent to converting an external mode to an internal mode (i.e. qqQ); 41 thus while fourth order scattering processes may still occur for external modes, any higher order processes are exceedingly unlikely. For the purpose of discussion in this work, we take .
Although the initial up-pumped vibrational energy remains localized on molecules, it will eventually dissipate spatially at a rate that depends on the anharmonic coupling between molecules. 42 In order for sufficient energy to be localized on an individual molecule, a large amount of energy must therefore up-pump at a rate that is faster than the dissipation rate. It follows that those materials which, according to Equation 4, exhibit faster vibrational up-pumping, can be expected to display greater susceptibility to initiation by mechanical impact.
The probability of scattering is proportional to the number of phonons available to scatter. Thus, the initial up-pumping events are dominated by qqQ processes, in which the inequality ( ) < applies. The states are known as the doorway ( ) + ( ) ( ) modes, as they act as efficient conduits for vibrational energy transfer between extended phonon modes and localised internal vibrational modes. This efficiency is largely a consequence of their mixed nature. As energy and phase must be conserved, the energy transfer step is dominated by processes in which . 29 ( ) = ( )/2 Once the doorway modes are excited, they can subsequently scatter to excite higher frequency modes in a 'ladder'-type approach. 43 The overall-process can be summarized by Figure 1.
In order to induce the initiation of a chemical reaction, a dissociation barrier must be overcome. In our theory this requires that sufficient vibrational energy is transferred from the lattice to localised molecular modes. The (sub)set of vibrational modes that are capable of inducing a bond dissociation are known as the target modes. The chemical simplicity of the binary energetic azides, M x (N 3 ) y , permitted identification of a single vibrational normal coordinate that is responsible for initiating . 37 In such cases, a N -3 direct up-pumping model can be discussed. 25 Such a mechanism is possible if either (1) the target mode is sufficiently isolated such that the multi-phonon density of states is negligible in its vicinity, or (2) the target mode has sufficiently strong coupling to the phonon bath to permit rapid, direct transfer. While condition (1) often applies to EMs composed of simple energetic molecules such as for other EMs, such as those consisting of cyclic molecules N -3 , shown in Figure 2, the decomposition pathways can be quite complex. Instead, an up-pumping model can be constructed that is based on an indirect (thermal) mechanism. 41 In this case, vibrational energy rapidly up-pumps into the entire manifold of internal molecular modes which leads to the population of a subset of vibrational modes responsible for bond rupture and initiation. This is the approach used throughout this work. Sensitivity is thus correlated with the total up-pumping into internal vibrational modes. This correlation is consistent with RRKM theory, in which localisation on particular modes becomes more probable with higher amounts of localised energy. Furthermore, the 0 K limit previously employed in our earlier work on the azide materials 37 is now extended to include temperature effects and includes the relative heating associated with adiabatic compression. Figure 1: Schematic representation of the vibrational up-pumping model. The lowest frequency modes (black rectangle) absorb the initial mechanical energy, which then transfers to the 'doorway' modes (red rectangle). Finally, the doorway mode energy uppumps into internal molecular modes at higher frequency to induce bond rupture. 1,1'-Azobistetrazole ABT ----< 1 44 4,5-Dihydro-5-nitrimino-1H-tetrazole DNIT ----1.5 45 Hexanitrobenzene HNB 11 2.75 -- 46 2-Methyl-5-nitramino-2H-tetrazole MNT ----3 45 1,3,5,7-Tetranitro-1,3,5,7-tetrazocane HMX 32 Recalculated for a drop weight mass of 2.5 kg. b Assumed value.

Test Set
Typical organic EMs consist of molecules with a large number of atoms and, common to other molecular materials, form crystals that adopt comparably large unit cells that possess low symmetry. Test candidates were therefore selected for which calculations were computationally tractable whilst ensuring that a broad range of impact sensitivities and structure types were represented. The test set was explicitly selected to include families of structurally related compounds, as well as a range of chemically diverse species. In this way, we aimed to test our model both within and between families of EMs, the latter having proved particularly challenging in the literature. Although many of these materials have been studied extensively, it is not uncommon for a range of impact sensitivities to be reported for the same compound. Documented factors influencing the reproducibility of impact sensitivity measurements for a given compound are the particle size, particle size distribution and phase purity. 52 Three of the compounds, guanidinium tetrazolate (GTZ), aminoguanidinium tetrazolate (AGTZ) and diaminoguanidinium tetrazolate (DGTZ) are new materials, reported here for the first time. The test set used in our study is given in Figure 2 (with crystal structures in the ESI S1) and Table  1, alongside a summary of reported impact sensitivities. The output of the computational model will be compared against the trend in experimental impact sensitivities. This trend in impact sensitivity is based on the following collated observations based on the open literature and new measurements (for abbreviations see Table 1 and Figure 2): 1. It is widely acknowledged that TATB is insensitive to impact. 12 Most reports state its impact sensitivity to be 'immeasurable'. Only a single source suggests that TATB has a sensitivity on the order of 50 J, 52 although no indication was provided as to its origin or the experimental method used. 2. NTO is generally regarded as a highly insensitive material. 52 3. Despite discrepancies in the size of the sensitivity value (ranging from ~20 J to 50 J), it is accepted that -FOX-7 is less sensitive than -HMX. 52 4. HNB is well-known to be highly sensitive and exceeds the sensitivity of -HMX. 50 5. The sensitivity of 5-substituted 1H tetrazoles is ranked according to R groups and follows the series CH 3 < NH 2 < N 2 HO 2 < Cl < NO 2 < N 3 < . 53 The Cl-substituted tetrazole N + 2 is reported to have an impact sensitivity > 40 J. It is reasonable to expect the NH 2 -substituted compound to be effectively insensitive. No data could be found to support the notion that 5-amino-1H-tetrazole exhibits any notable sensitivity. For the purpose of this study we have assumed a value of > 50 J. The corresponding crystal structures are given in the ESI S1.

Results and Discussion
Electronic Structure-Band Gap Criterion EM initiation requires rupture of a chemical bond, which itself requires population of an antibonding state. This has led to the socalled 'band gap criterion' for EM sensitivity prediction. 55 The HSE06 screened hybrid functional is a recognised standard for band gap calculations, although it is often too expensive for large molecular systems. Based on results obtained for a series of inorganic materials, the less demanding hybrid GGA functional, B3PW91, has instead been shown to predict fundamental band gaps ( ) in good agreement with experiment. 56,57 For the study of the EMs considered here, it is therefore of interest to compare HSE06 and B3PW91 band gap energies alongside a standard GGA functional such as PBE.
To the best of our knowledge, only limited experimental data regarding values are available for the materials studied here. Experimental UV-Vis spectra have been documented for -HMX. 58 While the fundamental band gap (i.e. the difference between the ionization potential and electron affinities) is in principle different to the optical band gap (which is stabilized by electron-hole interactions, and may be restricted by symmetry), the discrepancy is often small in solid state organic materials. 59 Hence, the UV-Vis spectra should offer a good indication when considering the validity of calculated values. Furthermore, it should be noted that both -FOX-7 and TATB are yellow powders, which indicates the presence of optical transitions in the region of about 2.7 eV.
The values calculated for the series of molecular energetic compounds are given in Table 2. The B3PW91 functional predicts values of that are consistently slightly higher than those obtained from HSE06, ranging from ( (HSE06)) +0.24 eV to ( (HSE06)) +0.34 eV. As expected, due to the absence of exact exchange in the PBE calculations, considerably lower values are returned in comparison with the higher-level functionals. The values for PBE-based calculations available in the literature (see Table 2) differ only slightly from those calculated here. This is most notable for -HMX, which may be attributed to other differences between the computational models such as the use of plane wave bases sets for the literature reports, which contrasts with the localised basis sets used in this work. In all cases, the G 0 W 0 calculations found in the literature suggest a band gap that is larger than that calculated by either B3PW91 or HSE06. Given the limited availability of experimental data, and taking into account the colour of TATB and FOX-7, it can be inferred that G 0 W 0 quasi-particle methods may overestimate the values of . This has been demonstrated previously for inorganic systems. 56 As has been noted for a set of energetic azide materials, 37 there is no noticeable correlation between band gap values and the reported impact sensitivity of the molecular EMs. This also holds for the materials tested here, see Table 2. Based on the B3PW91 or HSE06 calculations, the predicted sensitivity ordering would be NTO > -FOX-7 HNB TATB > ABT > DNIT > ATZ > -HMX > ≈ ≈ MNT > HBT > GTZ > AGTZ > DGTZ. A comparison with the experimental ordering deduced above shows clearly that this is incorrect. Similar to observations for the azide materials, 37 there is no evidence of any correlation between sensitivity and the nature of the band gap in a material (either direct or indirect band gap). The band gap criterion performs poorly even within classes of compounds. For example, within the series of tetrazole-based compounds, an ordering on the basis of should deliver ABT > DNIT > ATZ > MNT > HBT > GTZ > AGTZ > DGTZ. While this ordering does capture the general trend of sensitive vs. insensitive, the accuracy of ordering remains poor compared to the data in Table 1. Thus, in order to find a correlation with experimental impact sensitivities, investigations must probe beyond electronic band gap energies. 60 Table 2: Fundamental electronic band gaps ( / eV) in the crystalline molecular energetic materials, arranged in order of decreasing impact sensitivity, and labelled as direct (D) or indirect (I). All values calculated here are based on a localised basis set.
In order to ensure the validity of the computational method, the simulated vibrational spectra for a subset of molecular EMs were first compared to experimental inelastic neutron scattering spectra (INS): -HMX, -FOX-7, NTO and TATB. INS spectroscopy is particularly suited for this purpose, as (1) it is sensitive to the motion of light atoms (particularly hydrogen), (2) it is not restricted to quantum selection rules, and (3) spectra are straightforward to simulate from calculated eigenvectors and vibrational frequencies.
The simulated INS spectrum for one of the more sensitive test compounds, -HMX ( Figure 3) shows, in general, a good agreement with experiment. The frequencies in the low-energy region of the INS spectrum are well-reproduced. The simulation underestimates the top level of the phonon bath ( ) by only ca Ω 5 cm -1 . There appears to be a systematic underestimation of the vibrational frequencies in the cm -1 region of the ω > 200 spectrum, amounting to ca. 20 cm -1 . The intensities of the lowest frequency modes are not reproduced particularly well, which may be indicative of minor issues in generating the correct eigenvectors; we do not however expect this to make any significant difference to subsequent calculations in this work. A better agreement between simulation and experiment is achieved for -FOX-7, Figure 3. The simulated INS spectrum of NTO also exhibits excellent agreement with experiment, Figure 3, despite a small instability (negative frequencies) that could not be converged out of the calculated phonon dispersion curves (ESI S3). Both the intensities and frequencies appear well-reproduced across the spectrum. This strongly suggests that this minor instability has only a negligible influence on the vibrational structure. There appears to be some disagreement among the vibrational modes at higher frequency (ca. 800 cm -1 ) in the INS spectrum of TATB ( Figure 3); overall, however, the spectrum is reproduced very well in both frequencies and intensities. A description of the effects of q-vector sampling density on modelling the INS curves is given in ESI S4.

Vibrational Structure
The value of is defined as the top level of the manifold of Ω lattice-based vibrational modes (i.e. the external modes). This definition stems from the average anharmonic approximation, which assumes that all external modes (q) couple sufficiently quickly in order to permit rapid equilibration. Solid materials that contain molecules will have three acoustic bands and 6 -3 external modes. Large organic molecules typically exhibit lowfrequency internal vibrations at energies below those of the highest-frequency lattice mode, also known as amalgamated modes. Within the context of the model developed here, we therefore define Ω as the highest frequency band above the external modes, which 6 cannot be distinguished from the phonon bath by a gap in the phonon dispersion curves. Generally, this allows for the capture of modes with high anharmonicity, thus participating in the rapid equilibration of mechanical energy. The resulting values of Ω are given for each material in Table 3. Deducing the value of Ω for NTO is somewhat challenging. While the 6Z vibrational external modes fall below cm −1 , amalgamated modes = 170 arising from -NO 2 wagging motions are observed at 200 cm -1 . As it is not clear which value constitutes the top of the phonon bath, both are considered in the following treatment. Note that, as discussed previously, 37 only those vibrational bands are considered that correspond to the energetic molecule, while in the case of salts those which solely arise from any counter ions, are not. This follows on from knowledge of the scattering rate, as given by τ, Equation 5, or analogously, Equation 7.
Where refers to the two-phonon density of states. This Ω (2) assumption requires that uncorrelated vibrational modes (i.e. those residing on different molecules) do not typically scatter to any notable extent. 42 Due to the typically low level of dispersion in organic materials, we treat all phonons as Einstein phonons. Hence, all terms in Equation 7 are wave vector-independent. This allows an evaluation of Equation 7 (and related equations) that is based on the phonon density of states, as opposed to the phonon dispersion curves.
The data presented in Table 3 gives an initial qualitative insight into predicting the ordering of impact sensitivity, namely that an increase in the gap between and the first doorway mode Ω correlates with a decrease in impact sensitivity.

Doorway Mode Density
Mechanical impact leads to direct excitation of the lowestfrequency lattice modes. 25 With the exception of amalgamated phonon modes, these low-frequency vibrations exhibit external lattice character, and hence do not induce a distortion of the internal molecular structure. It follows that excitation of the internal modes is required to achieve vibrationally-induced transient metallisation (i.e. closure of the band gap) 64 and to populate anti-bonding states in the energetic molecule. This process triggers initiation.

Within the average anharmonic approximation in Equation 7
(3) depends on the anharmonic character of the scattering vibrational modes. Hence, the phonon bath (i.e. ) equilibrates (3) ∝ rapidly within picoseconds. 27 Following equilibration of the phonon states, vibrational up-conversion occurs via phononphonon coupling up to a maximum frequency of within the 2Ω first anharmonic approximation. This represents processes in which the relationship applies. Vibrational states at (3) ∝ frequencies are defined as doorway modes. Ω < ω < 2Ω Subsequent excitation beyond the first anharmonic approximation requires scattering between a doorway mode and a phonon mode, for which applies and which therefore occurs (3) ∝ somewhat more slowly.

The rate of phonon scattering in Equation 7 is proportional to
. This is the Bose-Einstein population of a given phonon at a ( , ) given temperature which can be used to predict the rate of energy transfer between the phonon bath and doorway modes via ( ) processes. These occur at the rate (Eq. 8) 29 ℎ→

Equation 8
ℎ→ ∝ ( + 1) ( + 1) [ ℎ ( ) -( )]; Eq. 8 has been simplified by some authors using the energy transfer parameter, 29 . Hence, describes the κ temperature-independent energy conversion rate between a phonon and a doorway mode. Since the ratio and the scattering rate Ω/θ τ at T = 0 are both approximately constant for organic materials, Equation 9 depends largely on the number of doorway modes, , only. 36 The probability of up-conversion into the doorway states therefore depends on the density of doorway states (i.e. the number of modes per unit cell volume), Figure 4. This probability is calculated as the density of doorway modes. It is normalised by the magnitude of the doorway region and 3N (see ESI S5). It is generally found that the more sensitive materials contain a higher density of vibrational states within the doorway region, whereas insensitive materials exhibit only minimal density in this part of the spectrum. Hence, despite yielding a rather poor ordering of relative impact sensitivities for the EMs investigated here, this qualitative model does appear capable of offering a broad classification (sensitive vs. insensitive) of each material. Its predictions agree well with previous studies that employed this model. 35,36 Figure 4: Correlating density of doorway states against experimental impact sensitivity, with the more sensitive materials clustered within the orange circle, and the less sensitive materials in the green circle. Arrows indicate that higher impact energies could not be measured, and so are approximate values only.

Overtone-Based Model
The term doorway mode describes modes that are accessible in the initial excitation from the phonon bath. Hence, for scattering processes , a doorway mode can be defined (3) ∝ ( ) generally as having a frequency . Hence Ω < ω < ( + 1)Ω the scattering event (the first overtone of order has overtone of order has ) is limited to . = 2 1 = 2 = 3 < 3Ω This differs from previous works which have assigned ca. 700 cm -1 as being the limiting frequency for initial excitation, despite considering overtone processes with . 32,35,36 The seemingly ≫ 2 arbitrary choice of 700 cm -1 as the upper limit for vibrational uppumping appears to have stemmed from experimental limitations of previous decades, rather than being a material-specific value. This issue has been noted previously. 35 To consider propagation of overtone up-pumping pathways the phonon density of states, is generated to the order, 37 (ω) subsequently normalized by . 36 A numerical example is ∫ ω (ω) provided in ESI S6 and S7. This is done in agreement with an indirect up-pumping mechanism, in which up-pumped energy equilibrates across the internal vibrational modes of the material. This also acts to remove size dependencies associated with normalization of the density of states to 3N.
Noting that as the rate of phonon-phonon scattering decreases with increasing order, we limit discussion to overtones of order , ≤ 2 i.e. the first two overtone pathways. The resulting trend is shown in Figure 5 against available experimental impact sensitivity. It is generally observed that the sensitive materials are characterized by large values of , while less sensitive materials exhibit increasingly smaller values in a nearly exponential fashion. Two materials are seemingly anomalous to this trend: DNIT and AGTZ. The former is predicted to be insensitive, while the latter is predicted to be sensitive (within a large uncertainty window). The opposite appears to be true in both cases.
Including overtone pathways of leads to poorer correlation > 2 against experimental impact sensitivity (ESI S8). This is in fact an encouraging observation, as a physical basis should not be found by including the low-probability higher-order scattering processes. It is important to note that this contrasts with earlier reports, which do include discussion of higher-order overtones. 32,35,36 However, the data sets explored in these earlier reports have simplified the discussion to Brillouin zone centre ( only) frequencies only and did not consider population conservation between resonant states. These factors may well account for the differences observed.

Two-Step Up-Pumping Model
Historically, predictive theories that are based on vibrational uppumping have varied between two base models. In the first, the overtone contributions are considered, 32,35,36 while the second model considers up-pumping based on a two-phonon density of states model. 30,34 The latter is primarily composed of combination pathways, in which . Our previous work has suggested that q 1 ≠ q 2 both processes may be important when considering a predictive model for multi-phonon up-pumping. 37 Within the remit of an indirect up-pumping model, we develop herein a two-layered approach within the first anharmonic approximation that combines both base models. In step (1) the first overtone ( is projected = 1) onto the doorway modes (i.e. to frequencies , Ω < ω < 2Ω according to 27

Equation 11
∝ Ω (2) in which is restricted to , and where and represent Ω (2) ω/2 ℎ the Bose-Einstein populations at temperature, , of the phonon and doorway frequency, respectively. We note that for the purpose of this model, consideration is restricted to the initial energy transfer and, hence, vibrational cooling is not considered. This has the effect of transferring excess phonon population onto the doorway modes.
Time-dependent calculations by Kim and Dlott 27 suggest that for a model molecular material, secondary up-pumping begins within picoseconds of the initial process of populating the doorway region. Hence, in step (2) of our up-pumping model we introduce a combination-based up-pumping process via 30

Equation 12
∝ Ω (2) in which the doorway population is now defined as the sum of the thermal population and that already up-pumped according to Equation 11. Thus, where the target mode also sits within the doorway region, the value of includes the initial up-pumped ω population. Restricting discussion to the first anharmonic approximation, the total up-pumped density is considered across the range and normalised by the PDOS within this Ω →3Ω region. This effects the normalisation of the up-pumped density by the number of states into which the density is pumped.
In accordance with Equations 11 and 12, this two-step model permits consideration of temperature. If the model is considered at an equilibrium temperature of 300 K, Figure 6 is obtained. There is generally excellent agreement between the predicted ordering and the experimental sensitivity, with two notable exceptions: DNIT and NTO. It is promising, however, to find that the model corrects the ordering of the sensitive compounds, with the highly sensitive materials (HNB, ABT and MNT) all predicted to be higher than -β HMX. HNB is again predicted to be considerably more sensitive than ABT, and NTO (with = 170 cm 1 ) is a clear outlier, Ω although the simulated up-pumped density falls dramatically when = 200 cm -1 ( Figure 6). While this higher value of does Ω Ω not rigorously adhere to the definition of stated above, the Ω potential error in vibrational frequencies (particularly in light of the imaginary calculated frequency) warrants its consideration. Figure 6: Correlation of the two-layered up-pumping model at 300 K against experimentally determined impact sensitivity. A dashed line is added to aid visualisation.
As has been discussed earlier, the increase in equilibrium temperature will depend on both compressibility and anharmonicity of the material through Equation 3. Upon instantaneous compression, the shock energy is deposited into the phonon bath and, as the phonon heat capacity, , is considerably ℎ lower than the total heat capacity, the quasi-equilibrium temperature achieved by the phonon bath, will greatly exceed θ ℎ , that achievable by the bulk solid. The phonon heat capacity is given by the Einstein model as roughly per molecule in the primitive 6 unit cell. Due to the amalgamation of some molecular modes into the phonon bath, the approximation that for a ℎ ≈ 6 + primitive cell of molecules and amalgamated modes holds. This therefore reflects the number of modes within the phonon bath, Table 4. It is reasonable to assume that the compressibility and anharmonicity of the molecular materials studied here will be similar in this respect. Therefore, the respective values of scale ℎ as , where is the total number of vibrational = / ℎ ⇒ = / modes below . Ω As a roughly representative value for shock temperature a value of = 3000 K was chosen. When Figure 6 is re-examined with the θ ℎ addition of a shock temperature, the correlation with experimental impact sensitivity improves substantially, see Figure 7. At this stage, the sensitive EMs are, for the first time, correctly ordered, with a steep descent towards the insensitive materials. Importantly, NTO is placed correctly, but only if the higher-lying -NO 2 wagging vibrational modes are included in the phonon bath, i.e. when the higher value of is used. Furthermore, it becomes possible to Ω predict an ordering for the tetrazolate salts, whose experimental sensitivities have not been measured above 50 J.
It is worth highlighting that the correlation based on a multi-phonon up-pumping model given here is the strongest one established to date. Moreover, this model appears capable of ranking EMs across a chemically diverse set of energetic molecules. This finding suggests that all of the information regarding bond strength (i.e. dissociation energies) is captured within the relative vibrational structures of these materials and can thus be taken account of in a multi-phonon model. While many of the data points in Figure 7 appear to fit a decay line, it should be noted that, within the model proposed here, a completely insensitive material must display an up-pumping density of zero. Hence, it is reasonable to expect an asymptotic function to fit the data points and this has been explored in Figures 5 through 7. The validity of this type of correlation strengthens the suggestion that the value found Ω = 200 cm -1 for NTO is consistent across all datasets in this work. This finding is consistent with earlier works on HMX, which assume amalgamation of low-lying NO 2 waging modes into the phonon bath. 28 Figure 7: Up-pumping densities according to the two-layered multi-phonon approach for predicting impact sensitivity. The value of is set according to Table 4, with all modes set at θ ℎ ω > Ω . The dashed line is added as a visual aid. Arrows = 300 K indicate uncertainty in experimental sensitivity measurement based on the experimentally determined lower limits.

Towards Molecular Level Rationalisation of Impact Sensitivity Behaviour
The ability to up-pump vibrational energy depends critically on the density of vibrational states within the doorway region, Ω < ω . In the cases studied here, these vibrational modes < 2Ω manifest typically as ring deformation (for cyclic EMs) and the wagging motion of 'dangling' moieties. These doorway modes are soft internal molecular modes. It follows that molecules bound by weaker interatomic potentials (e.g. in the case of highly electron deficient systems) will contain a larger number of these doorway states. However, due to the intimate relationship between the crystal packing (via the phonon bath and ) and both the Ω position and size of the doorway region, it is not yet obvious how doorway density relates directly to molecular structure. At present, any correlation can only be understood in systems which share crystalline packing structures, such as TATB and HNB, Figure 8. It is a commonly assumed observation that layered materials are generally insensitive to mechanical impact. HNB is however known experimentally to be substantially more sensitive to impact than TATB. Interestingly, all of the up-pumping models presented in this work correctly predict the impact sensitivity ordering of these two compounds, thereby presenting an opportunity to determine why two structurally similar compounds exhibit such different impact sensitivity behaviour. Firstly, we note that the doorway density of HNB is considerably higher than that of TATB (Figure 4), which suggests that the former is less rigid than the latter. This region is dominated by the -NO 2 waging motions and ring deformation modes. The same motions in TATB occur at higher frequencies as the -NO 2 groups are held in the plane of the ring by internal hydrogen bond interactions with adjacent NH 2 groups. They are therefore outside this crucial vibrational region (see ESI S9). Secondly, we observe that the localised vibrational modes for HNB appear at lower wavenumbers than for TATB (see ESI S9), such that more modes fall below the upper limit for 3Ω vibrational up-pumping for the more sensitive compound. This again suggests that HNB is less rigid than TATB. The weakening in ring vibrational modes also correlates with a weakening in ring aromaticity in HNB when comparing the electron density distribution with that for TATB ( Figure 9). Thus, the vibrational up-pumping model presents the first indication of a physical mechanism that is able to rationalise the enormous difference in the sensitivity of two layered EMs. This finding has two important implications for the rationalising the design of new insensitive EMs. Firstly, in order to reduce the density of states in the doorway region and shift the localised vibrational modes that receive the up-pumped vibrational energy, rigid, electron dense molecules are required. This is characteristic of molecules such as FOX-7, NTO and TATB. Secondly, the value of should be minimised in order to minimise the upper bound Ω up-pumping limit of , as defined by the second overtone 3Ω population limit. This can be achieved by increasing the degree of layering in the crystallographic packing. The present contribution makes clear that the value of is not the crucial design Ω parameter, but it can assist to reduce sensitivity, provided that the density of the doorway states remains low.

Conclusions
In the present work, a vibrational up-pumping model is developed to rationalise the sensitivity of a chemically and crystallographically diverse range of molecular energetic materials (EMs). The fully ab initio vibrational structures for a subset of these EMs are in excellent agreement with experimental inelastic neutron scattering spectra (INS). This confirms the validity of the theoretically predicted vibrational structures. A qualitative analysis of the vibrational doorway mode density appears to offer an excellent means for ready assignment of sensitivity. However, these assignments are, at best, limited to 'sensitive' vs. 'insensitive'. Subsequently, vibrational energy up-pumping was explored using an overtone mechanism, where is restricted to the coupling ( ) of two counterpropagating phonons. As suggested in previous works, 36 this type of model does offer good agreement with experimental sensitivities, despite failing in some cases. However, the correlation can be improved considerably when a two-layer up-pumping mechanism is adopted. This mechanism implies an initial excitation of the doorway states that occurs by overtone pathways. Subsequently, the overtone pathways are pumped higher via combination pathways with the underlying vibrational density of states. This produces excellent agreement between predicted up-pumping rates and experimental sensitivities. This is particularly striking when differences in the relative shock (adiabatic heating) temperatures are considered. Despite its ability to predict the relative impact sensitivity of a broad range of EMs, it must be remembered that the model remains limited to ideal crystal structures. The model predicts the intrinsic sensitivity of EMs, but cannot yet describe the influence of crystallographic defects, grain size, distribution and grain boundary effects. Nevertheless, the new model is still of considerable benefit for screening new materials, and indicative of their inherent potential towards initiation, once the crystal structure is known.
Furthermore, it is important to highlight that the work presented here is based on the internal flow of energy which follows an external perturbation. For the purposes of this work, we have assumed that the excitation occurs via a conserved adiabatic heating process. Alternative excitation processes, including fracture and contact stresses (e.g. with grit) have also been suggested. Importantly, the nature of the excitation primarily affects the magnitude and spatial distribution of the excitation. The model presented here is largely independent of these parameters. Further work is instead required to understand the magnitude of excitation imposed by these alternative mechanisms and will be the subject of follow-up investigation.
As a result of this work it is now possible to suggest design targets for new EMs. Foremost in the search for insensitive materials appears to be the need for a low density of doorway vibrational modes and for the localised molecular vibrations to be shifted to higher wavenumbers. Both of these features, which are molecular design targets suggest that increasing structural rigidity will suppress vibrational up-pumping. An excellent example of this is provided in the comparison of the structurally very similar compounds TATB and HNB. The value of also appears to be crucial for suppressing rapid uppumping into the highest levels of the vibrational structure.
This design parameter falls within the realm of crystal engineering, and suggests a rationalisation for the differing sensitivities of polymorphic and multi-component crystal forms encountered previously.
The ab initio model developed here offers excellent predictive capacity across a chemically and crystallographically diverse range of EMs. By introducing consideration of temperature and heat capacity, it allows further exploration of temperature dependence, and dynamic effects associated with material compressibility upon impact. Further development of this model therefore offers immense potential to rationalise EM property behaviour in terms of experimentally tuneable parameters.

Experimental and Computational Methods Preparative Procedures and Analytical Chemistry
NTO. 3-nitro-1,2,3-triazol-5-one (NTO) was prepared by nitration of TO (prepared by Dr. S. Kennedy, School of Chemistry, University of Edinburgh), according to Scheme 4.1. 65 TO (ca 5 g, 60 mmol) was dissolved in excess concentrated nitric acid (70%), and held at 55 o C for ca 45 minutes. The nitration was quenched in an ice bath, filtered and rinsed in ice water. Purity was verified by X-ray powder diffraction.
TATB and -HMX. Samples of these materials were taken from β available stock. -HMX was provided by the Cavendish β Laboratories, University of Cambridge. TATB was synthesised by C. Henderson (School of Chemistry, University of Edinburgh). Samples were used as provided, without further purification.
The crystallographic information files for compounds GTZ (CCDC 1900007), AGTZ (CCDC 1900008) and DGTZ (CCDC 1900009) are available free of charge from the Cambridge Crystallographic Data Centre (CCDC). 51 A structure discussion of the crystals of the compounds GTZ, AGTZ, DGTZ (reported here for the first time) as well as triaminoguanidinium tetrazolate (TGTZ) will be provided in a forthcoming paper.

Condensed Matter Calculations.
The PBE-D2 scheme has been previously demonstrated to work well for molecular energetic materials. 66 Optimisation was performed using plane wave DFT as implemented in CASTEP v16. 74 The electronic wavefunction was expanded in plane waves to a kinetic energy of 1800 eV for all systems, except HBT (cut-off 1600 eV), and -HMX and TATB β (cut-off 1300 eV). Forces were converged to < eV Å -1 , 5 × 10 -4 and stresses to < GPa. The energy change per atom was 5 × 10 -4 converged to < eV/atom. The resulting unit cell 1 × 10 -9 parameters are given in Table S3.1. In all cases, norm-conserving pseudopotentials were used, as available within the CASTEP v16 software package. All phonon calculations were based on primitive unit cells. Phonon calculations were performed using the linear response method to calculate the dynamical matrices on a regular grid of wave vectors, and Fourier interpolated to a fine grid of > 150 points for generation of phonon density of states (DOS, )). (ω Phonon dispersion curves were generated by Fourier interpolation of the computed dynamical matrices along high symmetry paths as proposed by SeeKPath 75 and labelled by IUPAC convention. Electronic band structures were calculated in CRYSTAL17 76 using localised basis sets, available from the CRYSTAL17 database. Their selection was based on previous success with similar materials and DFT functionals (H-H_pob_TZVP_2012 77 ; C-C_m-6-311G(d)_Heyd_2005 78 ; N-N_m-6-311G(d)_Heyd_2005 78 ; O-O_m-6-311G(2d)_Heyd_2005 78 . In order to ensure closest reproduction of experimental results, all calculations were performed on the crystallographically determined geometries. The band structures were calculated using the HSE06, 78 B3PW91 79 and PBE 80 functionals. For all materials, the tolerances (TOLINTEG) were set at 7 7 7 9 30 (as recommended for use with these basis sets 78 ). The electronic structure was sampled across a regular grid of points, with ca 120 points sampled in each material.
Inelastic Neutron Scattering Spectroscopy. All INS spectra were collected using the TOSCA spectrometer at the ISIS Neutron and Muon source. 81,82 Samples (ca 1.5 g) were placed in aluminium sample holders. Samples were cooled to ca. 10 K and data collected for a total of ca. 400 Ah. The sample temperature was μ subsequently increased in steps of 50 K up to a maximum of 200 K for -HMX and TATB, and 150 K for -FOX-7. Data were β α collected at each 50 K interval. Both forward and back-scattered data were accumulated and corrected for scattering from the sample holder and background. All data processing was done using Mantid. 83 INS spectra were simulated using ABINS, 84 as implemented in Mantid. Only first-order quantum events (i.e. the fundamentals) are considered in the simulation of INS spectra.