The development of new energetic materials (EMs) is accompanied by significant hazards, prompting interest in their computational design. Before reliable *in silico* design strategies can be realized, however, approaches to understand and predict EM response to mechanical impact must be developed. We present here a fully *ab initio* model based on phonon up-pumping that successfully ranks the relative impact sensitivity of a series of organic EMs. The methodology depends only on the crystallographic unit cell and Brillouin zone center vibrational frequencies. We, therefore, expect this approach to become an integral tool in the large-scale screening of potential EMs.

## INTRODUCTION

The mechanically induced reactivity of solids represents a fascinating and poorly understood aspect of materials science. With respect to energetic materials (explosives, propellants, pyrotechnics, and gas generators; EMs), this reactivity is crucial for determining their handling safety and potential for real-world applications in both civilian and military technologies.^{1} EMs are characterized by their ability to rapidly convert large amounts of chemical potential energy into kinetic energy, with potentially devastating effects if uncontrolled. Owing to stringent safety and performance requirements (e.g., detonation pressure and velocity and the heat of explosion), only a limited number of EMs have found widespread applications. However, many of the well-established EMs do not comply with increasing environmental and public health regulations, prompting the need to develop new, non-toxic, and “green” EMs that compromise neither on safety nor on performance.^{2}

The impact sensitivity (IS) of an EM represents an internationally recognized material hazard classification criterion.^{3} Where materials are highly sensitive (i.e., initiate on mild impact), they are generally not fit for real-world applications, regardless of their other performance characteristics. Correspondingly, significant efforts have been devoted to identifying strategies for characterizing, rationalizing, and designing low sensitivity EMs (LSEMs) that offer greater potential for safe handling. Various experimental protocols exist to measure the IS of EMs.^{4} These approaches typically involve the use of a drop hammer device, such as the BAM Fall-hammer or the Rotter impact device, both requiring the preparation of bulk quantities of novel EMs. Thus, while the discovery of new EMs requires significant investment of both time and financial resources, the lack of *a priori* knowledge regarding their explosive properties also carries considerable risk to health and safety.

The challenges presented by the experimental discovery of LSEMs have prompted efforts to derive theoretical descriptions of EM sensitivity.^{5,6} Many efforts have focused on the screening of molecular parameters. For example, bond dissociation energies have been successfully correlated against EM sensitivity,^{7} including through machine learning analyses.^{8,9} A correlation between the IS and molecular oxidizing potential (i.e., oxygen balance) has been successfully demonstrated for a series of polynitroaromatic explosives,^{10} although attempts to verify the same correlation for a series of PETN analogs failed.^{11} Other molecule-centric approaches have included the correlation of sensitivity with NMR chemical shifts,^{12} atomic composition and their corresponding HOMO–LUMO gaps,^{13} electronegativity,^{14} and charge density distribution.^{15}

While these approaches have proved effective for rationalizing sensitivity trends in chemically related series of EMs, many important phenomena remain outside their predictive capacity. Reports have shown, for instance, that IS can vary for the same EM expressed in different polymorphic forms. For example, HMX (octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine, see Fig. 1) exists in the thermodynamically stable *β*-form under ambient conditions, but it is often contaminated with traces of the α- or γ-forms, which are both more sensitive to mechanical impact.^{16,17} An even more sensitive form (δ-HMX) is readily accessible upon heating.^{18} Similarly, a growing number of studies have shown that the formation of multi-component crystals (salts and co-crystals) can have significant influence on IS properties.^{19,20} The crystal lattice, therefore, must play a role in determining IS, and consequently, crystal-centric approaches have become the focus of growing attention. To this end, IS correlations have been sought against a number of solid-state characteristics, including intermolecular bond distances,^{21} electronic band gaps,^{22} crystal void space,^{23,24} and crystal morphology.^{25} However, these approaches generally lack physical insight into the initiation mechanism and have, therefore, been largely engulfed by more sophisticated models. For example, identifying the need to associate bulk compressibility with electronic excitation, Bondarchuk successfully correlated IS against pressure-induced metallization for metal azides^{26} and across a breadth of molecular EMs.^{27} In an alternative, albeit related, approach, various authors considered the development of local stress in perturbed solids, which results from anisotropic distortion of the solid under mechanical loading.^{28,29} Following on from these efforts, correlations between crystal packing geometry (e.g., *π*⋯*π* stacking)^{30,31} have become increasingly popular.^{28,32}

Following the pioneering work of Dlott and Fayer,^{33} we recently presented an *ab initio* phonon up-pumping based approach for predicting EM IS.^{34,35} This up-pumping model captures critical features of the impact event, including compression-induced excitation^{36} and vibrationally induced metallization.^{37} Our previous model was based on computationally expensive phonon dispersion curves (PDCs). The significant time and resources required to obtain PDCs represent a substantial barrier to implementing our model for larger-scale material screening. Moreover, the need for PDCs places large energetic molecules in low symmetry unit cells out of reach. A recent study by Bernstein^{38} suggested that up-pumping based approaches using limited sampling of the Brillouin zone can be sufficient to achieve positive correlation. This has motivated our current study to explore the potential to expedite our computational screening approach by considering only the Brillouin zone-center (Γ-point) *ab initio* vibrational frequencies. A key consideration in our up-pumping model is the assignment of Ω_{max}, the vibrational frequency denoting the change in the mode character from lattice to molecular behavior. In this work, we explore assignment through three separate methods, namely, mode counting, consideration of the vibrational mode contribution to the heat capacity, and finally, through tracking the changes in the molecular center of mass expressed by the eigenvectors. Herein, we briefly outline the theoretical background of our phonon up-pumping model and apply it to a well-known structurally and energetically diverse set of EMs.

## THEORY

The mechanical compression of a material leads to a bulk increase in the thermodynamic energy Δ*U*. This energy distributes between heat *H* and work *W* such that Δ*U* = *H* + *W*. Partition of the input mechanical energy, and hence, the elevation of heat in the system, can be approximated from consideration of the compression factor *V*/*V*_{0} and the Grüneisen parameter, γ. The final temperature *T*_{f} that is reached by an adiabatically compressed solid at initial temperature *T*_{0} with initial volume *V*_{0} can be described according to^{40}

The total heat transferred to the system can be re-written in terms of the bulk material heat capacity, *C*(*T*), with

Upon mechanical impact, a compressive wave passes through the material, coupling most strongly to the lattice acoustic modes. Impact energy, therefore, transfers into these few branches.^{41} Owing to the significant phonon–phonon scattering that occurs between lattice phonon modes, this energy rapidly redistributes throughout the external acoustic and optical branches. A good approximation for the initial conditions directly following mechanical impact is, therefore, to assume that the energy from the impact event is distributed across all of the lattice modes, together referred to as the phonon bath.^{33} The initial phonon bath quasi-temperature *ϕ*_{ph} condition can then be described by replacing *C*(*T*) with the phonon heat capacity *C*_{ph} such that

assuming that the phonon quasi-temperature must be a state of the bulk temperature–volume curve.^{33} As the phonon bath represents only a small number of the total vibrational states in the material, the phonon heat capacity is considerably lower than the bulk heat capacity. Correspondingly, the phonon quasi-temperatures reach significantly higher values than will be achieved after thermalization of the material. Hence, the initial conditions for our model comprise a vibrationally “hot” phonon region and a vibrationally “cold” molecular region of vibrational states.

An EM initiation event results from the rupture of a covalent bond. The energy stored in the phonon bands must, therefore, up-convert to reach an active molecular vibrational state, which is termed a target mode. To this end, we consider the vibrational Hamiltonian of a molecular crystal,

where *H*_{i} describes the internal motion of molecules, *H*_{L} denotes the motion of molecules about their equilibrium positions (i.e., the lattice vibrations), and *H*_{c} provides a route for coupling between the internal and external modes. For the general case of indirect phonon up-pumping, *H*_{c} can be expanded as

where

For simplicity, we have assumed Einstein phonons in Eq. (5), thereby avoiding explicit consideration of wave vector dispersion. This is generally a good approximation for the internal vibrations of molecular materials. The Kronecker delta of Eq. (5), which ensures momentum and energy conservation, is conveniently rewritten as *ρ*^{(2)}. This is the two-phonon density of states and describes the formation of a new phonon state Ω_{i} from the scattering of two lower energy phonons *ω*_{j} and *ω*_{k}. The coefficient *V*^{(3)} describes the anharmonic coupling strength between phonons *i*, *j*, and *k*. Spectroscopic evidence^{42} and semi-empirical calculations^{43} have suggested that the values of *V*^{(3)} are approximately constant for molecular energetic materials. Correspondingly, the relative scattering depends only on the function *ρ*^{(2)}. Early developments^{33} demonstrated that vibrational up-pumping follows a multi-staged approach. Energy from the “hot” phonon bath (i.e., from *ω* < Ω_{max}, where Ω_{max} denotes the top of the phonon bath) is first transferred to the doorway states (denoted by Ω_{max} < *ω* < 2Ω_{max}) with a rate given by *τ*_{1},

Here, it is assumed that *τ*_{1} transfers all excess population energy (denoted by the Bose–Einstein populations, *n*) from the phonon states (*n*_{ph}) into the doorway region (*n*_{d}), thereby re-establishing equilibrium populations of *n*_{ph} prior to further up-pumping. The “hot” *n*_{d} states subsequently up-pump to the target modes (*n*_{T}) in the second step according to rate *τ*_{2}. Where dispersion of *ω*(** k**) is large (i.e., in the phonon bath), the Einstein approximation [

*ω*(

**) =**

*k**ω*] is not suitable. This leads to the energy restriction in Eq. (6a). Where the target mode (

*n*

_{T}) is adequately described as an Einstein phonon, explicit consideration of wave vector conservation can be neglected. This is done in Eq. (6b). For simplicity, we can borrow the terms

*overtone*and

*combination*from vibrational spectroscopy to describe the analogous processes involved in Eqs. (6a) and (6b), respectively. In this terminology, we note that the contributions to

*ρ*

^{(2)}in Eq. (6a) are limited to the first overtone, as described by

*ω*=

*ω*

_{T/2}.

All up-pumping steps are assumed to require delocalized vibrations, i.e., to include at least one phonon bath mode. Correspondingly, we generate *ρ*^{(2)} in Eq. (6) by assuming that at least one phonon must fall at *ω* < Ω_{max}. Similarly, the second phonon mode must not be in thermal equilibrium with the bulk and must, hence, fall at *ω* < 2Ω_{max} in Eq. (6b). The highest available excitation by our phonon up-pumping model is, therefore, 3Ω_{max}.

The up-pumping model developed here is reminiscent of the starvation kinetics model developed by Eyring for unimolecular reactions at a shock front.^{39} In the Eyring model, a dissociating bond (i.e., the target bond) is at equilibrium with a “cold” (or “starved”) reservoir of oscillators (analogous to doorway modes in our model). As the shock front propagates through the material, it does not equilibrate directly with the dissociating bond but rather transfers energy to the starved reservoir, which, in turn, “feeds” the dissociating bond (*cf.* our two-step up-pumping mechanism). Although the conceptual construction of our model is consistent with the Eyring model, the formalisms differ. Most notably, the Eyring model derives its system dependence by considering bond dissociation energies, thereby requiring the knowledge of the primary decomposition event. In contrast, system dependence in our model follows from the relative vibrational energy transfer pathways, without any need for *a priori* knowledge of the material reactivity. We envision, however, the future need to include into our model concepts developed by Eyring.

## COMPUTATIONAL DETAILS

### Gas phase simulations

All gas phase simulations were performed in ORCA v4.2.^{44} The molecular structures were extracted from available crystallographic data and fully relaxed at the PBE/def2-TZVP level. The force constant matrix was calculated analytically for each relaxed geometry. Owing to the known deficiency of this level of theory at reproducing experimental frequencies, we adopt the corresponding scaling constant 1.0306.^{45}

### Simulations of periodic systems

All periodic DFT calculations were performed within the CASTEP suite v 19.11.^{46} All input structures were taken from the Cambridge Crystallography Data Centre (CCDC) database: Hexanitrobenzene (REF: HNOBEN), *ϵ*-CL-20 (REF:PUBMUU12), *β*-HMX (REF: OCHTET13), m-TNT (REF: ZZMUC), *α*-FOX-7 (REF: SEDTUQ01), *α*-Nitrotriazolone (REF: QOYJOD06), and TATB (REF: TATNBZ03). In each case, the electronic structure was expanded in plane waves to a maximum cut-off energy of 1200 eV and the nuclear–electron interactions were approximated with norm-conserving pseudopotentials as generated on-the-fly in CASTEP. The Brillouin zone was sampled on a Monkhorst–Pack grid^{47} with spacings no greater than 0.05 Å^{−1}. Owing to its proven success on these materials,^{35,48,49} the exchange-correlation functional of Perdew–Burke–Ernzerhof (PBE)^{50} was used in all cases and the Grimme D2 dispersion correction.^{51} The structure was considered fully relaxed when the convergence of residual forces <1 × 10^{−4} eV Å^{−1} and the convergence of the SCF cycles <1 × 10^{−13} eV atom^{−1} were achieved. All phonons were calculated at the Brillouin zone center (**k** = 0) according to linear response theory, as implemented in the CASTEP suite.^{52} The acoustic sum rule was imposed analytically. LO-TO splitting was not considered.

## TEST SYSTEMS

In our previous work, an up-pumping model based on sampling of the complete phonon dispersion curves was presented.^{35} For comparison, herein, we select a subset of the same systems, namely, HNB (hexanitrobenzene), *β*-HMX (1,3,5,7-tetranitro-1,3,5,7-tetrazoctane), *α*-FOX-7 (1,1-dinitro-2,2-diaminoethene), *α*-NTO (nitrotriazolone), and TATB (1,3,5-triamino-2,4,6-trinitrobenzene). Additionally, our test set includes two new systems that crystallize in large unit cells, *ϵ*-CL-20 (hexanitrohexaazaisowurtzitane; REF: PUBMUU12) and the orthorhombic polymorph of TNT (2,4,6-trinitrotoluene; REF: ZZZMUC01), o-TNT.

The selection of molecules is chosen to cover a range of impact sensitivities, from the highly sensitive HNB through to the insensitive TATB. Moreover, the selected molecules represent a range of so-called “trigger linkages” (i.e., active explosophoric bonds), thereby assessing the generality of the approach. The sensitivity values for each compound were taken from the literature, Table I. Our model is based on the consideration of the crystallographic structure, with key parameters of our test system provided in Table II. We note that reports of impact sensitivity values vary widely in the literature, owing to differences in experimental conditions, including the particle size^{53} and testing conditions. Where possible, a range of sensitivity values has been provided, although we note the relative ordering of compounds does not depend significantly on this variability. We note that the sensitivity of o-TNT has not been widely reported. However, with near identical crystallographic and molecular structures, density, and compressibility as compared with m-TNT,^{54} we take its IS value to be equivalent.

EM . | IS (Nm) . | V_{exp} (Å^{3})
. | V_{calc} (Å^{3})
. | ΔV (%)
. | Space group . | References . |
---|---|---|---|---|---|---|

HNB | 2.75 | 581.49 | 585.25 | +0.6 | I2/c | 55 |

ϵ-CL-20 | 3–4.25 | 1397.25 | 1461.01 | +4.6 | P2_{1}/c | 56 |

β-HMX | 6.5–8 | 501.37 | 516.01 | +2.9 | P2_{1}/c | 56 and 57 |

o-TNT | 24 | 1823.55 | 1839.64 | +0.8 | Pca2_{1} | 58 |

α-FOX-7 | 25.2–31.8 | 522.33 | 514.32 | −1.5 | P2_{1}/n | 58 and 59 |

α-NTO | 72.75 | 902.06 | 927.22 | +2.8 | $P1\xaf$ | 57 |

TATB | 122.5 | 425.25 | 434.69 | +2.2 | $P1\xaf$ | 57 |

EM . | IS (Nm) . | V_{exp} (Å^{3})
. | V_{calc} (Å^{3})
. | ΔV (%)
. | Space group . | References . |
---|---|---|---|---|---|---|

HNB | 2.75 | 581.49 | 585.25 | +0.6 | I2/c | 55 |

ϵ-CL-20 | 3–4.25 | 1397.25 | 1461.01 | +4.6 | P2_{1}/c | 56 |

β-HMX | 6.5–8 | 501.37 | 516.01 | +2.9 | P2_{1}/c | 56 and 57 |

o-TNT | 24 | 1823.55 | 1839.64 | +0.8 | Pca2_{1} | 58 |

α-FOX-7 | 25.2–31.8 | 522.33 | 514.32 | −1.5 | P2_{1}/n | 58 and 59 |

α-NTO | 72.75 | 902.06 | 927.22 | +2.8 | $P1\xaf$ | 57 |

TATB | 122.5 | 425.25 | 434.69 | +2.2 | $P1\xaf$ | 57 |

. | . | . | . | Ω_{max}/cm^{−1}
. | |||
---|---|---|---|---|---|---|---|

EM . | Z/unit cell
. | N/molecule
. | Y
. | MC . | EV-MC . | CoM . | HC^{a}
. |

HNB | 2 | 24 | 9 | 215 | 140 | 140 | 145 |

ϵ-CL-20 | 4 | 36 | 16 | 217 | 217 | 217 | 180/222 |

β-HMX | 2 | 28 | 11 | 237 | 191 | 191 | 106/196 |

o-TNT | 8 | 21 | 9 | 212 | 212 | 212 | 160/217 |

α-FOX-7 | 4 | 14 | 3 | 169 | 169 | 169 | 174 |

α-NTO | 8 | 11 | 2 | 251 | 199 | 199 | 166/204 |

TATB | 2 | 24 | 6 | 140 | 140 | 140 | 145 |

. | . | . | . | Ω_{max}/cm^{−1}
. | |||
---|---|---|---|---|---|---|---|

EM . | Z/unit cell
. | N/molecule
. | Y
. | MC . | EV-MC . | CoM . | HC^{a}
. |

HNB | 2 | 24 | 9 | 215 | 140 | 140 | 145 |

ϵ-CL-20 | 4 | 36 | 16 | 217 | 217 | 217 | 180/222 |

β-HMX | 2 | 28 | 11 | 237 | 191 | 191 | 106/196 |

o-TNT | 8 | 21 | 9 | 212 | 212 | 212 | 160/217 |

α-FOX-7 | 4 | 14 | 3 | 169 | 169 | 169 | 174 |

α-NTO | 8 | 11 | 2 | 251 | 199 | 199 | 166/204 |

TATB | 2 | 24 | 6 | 140 | 140 | 140 | 145 |

^{a}

Multiple values reported due to the onset of plateau regions in the HC model below the MC-Ω_{max} value.

## RESULTS AND DISCUSSION

### Vibrational spectra of energetic materials

Prior to considering a vibrational up-pumping model, it is necessary to identify the upper limit of the phonon bath, Ω_{max}. This value is intended to bound all lattice modes and must refer to a point where *g*(*ω*) = 0 in the phonon density of states (see Fig. 2). While direct visualization of the eigenvectors (EV), in principle, offers a straightforward approach to marking the distinction between the lattice mode and molecular mode behavior, in practice, ambiguities often prevail.^{35} In molecular crystals comprising *Z*, *N*-atom molecules, there are six degrees of freedom (translations, rotations) for each molecule in the unit cell. It follows that, each wave vector in a molecular crystal comprises *6Z* external modes of vibration. Within the rigid body approximation, the *6Z* external modes are assumed to be well separated from the remaining *3NZ-6Z* internal molecular modes. However, a small number of internal modes—typically rocking and wagging motions of, e.g., NO_{2} moieties—can often appear at frequencies that are on par with the external modes. These are known as amalgamated modes. The total number of vibrational bands that are bound below Ω_{max} is, therefore, *Z*(6 + *Y*), where *Y* is the total number of amalgamated modes per molecule.

The first approach we adopted to identify an appropriate value for Ω_{max} was mode-counting (MC), wherein a direct comparison between the isolated molecule vibrational spectrum and its corresponding condensed phase density of states was made. The lowest frequency *6Z* external modes can be readily identified in the condensed phase calculation, thereby defining a maximum cut-off frequency for the external mode frequencies, Ω_{ext}. From our gas phase simulations, the number of internal modes that fall below Ω_{ext} are subsequently identified, *N*_{int}. Owing to the symmetry splitting of molecular modes in the solid state, this yields *Y* = *ZN*_{int} amalgamated modes. As these *Y* amalgamated modes do not contribute to the required *6Z* external modes, we extend the limit of Ω_{ext} to include the additional *Y* vibrational frequencies. Additional molecular modes that fall below the extended Ω_{ext} are included, and the value of Ω_{ext} is again updated. This procedure is repeated until no further molecular modes can be identified below the *6N* external mode maximum limit. The final cutoff is taken to be the top of the phonon bath, Ω_{max}. Values of Ω_{max} identified using this MC procedure are tabulated in Table II and discussed further in the supplementary material. It is worth highlighting that the values of Ω_{max} obtained in this way are generally consistent with our previous report based on direct visualization of the EVs.^{35} In a few cases (notably, HNB, β-HMX, and α-FOX-7), crystal field effects lead to marked shifts in the gas phase frequencies. Where this is the case, EV visualization can then be invoked to facilitate the tuning of MC-derived Ω_{max} values (see Table II and the supplementary material). In practice, the absolute eigenvalues obtained for both MC- and EV-MC-Ω_{max} are extended by 5 cm^{−1} when applied to the Gaussian-broadened phonon density of states.

To facilitate EV analysis, we note that internal modes should be characterized by negligible displacement of the molecular center of mass (CoM), which contrasts with the marked displacement of the CoM expected for external modes. Correspondingly, we explored the possibility for rapid EV analysis—and hence, identification of Ω_{max}—by tracking the variation in the CoM as expressed by the 3*N* rectilinear eigenvectors. In this way, we obtain Ω_{max} values as the absolute vibrational eigenvalues above which no eigenvectors display CoM displacements (see the supplementary material). In practice, 5 cm^{−1} are again added to these values to allow for Gaussian broadening. The results obtained are essentially identical to the EV-MC method, which was as expected as both involve direct analysis of the eigenvectors. This route does, however, provide a quantitative analysis that is much less subjective than the visual comparison inherent in the EV-MC method and also provides a route to the automated screening for Ω_{max} values.

With the limiting value Ω_{max} set by the EV-MC (or analogously CoM) method, it becomes straightforward to calculate the heat capacity values required to solve Eq. (3). Based on the Brillouin zone center vibrational frequencies only, an approximate density of states was generated by applying a Gaussian smearing of ±5 cm^{−1}. This yields a FWHM for isolated vibrational bands of ±5 cm^{−1}.^{42} We, therefore, obtain the bulk heat capacity by Eq. (3) and the phonon bath heat capacity via the same equation by imposing an upper integration limit of Ω_{max} (see the supplementary material). It follows from Eq. (3) that the adiabatic heating in each case becomes proportional to the ratio *C*_{tot}/*C*_{ph}. The high temperature limit of the bulk and phonon heat capacities can then be determined according to

where *ω*_{i} = Ω_{max} for the phonon heat capacity and *ω*_{i} extends across all *3N* vibrations for the total heat capacity. Based solely on the zone center vibrational frequencies, it is worth noting that the heat capacities are expected to be over-representative for each system (see the supplementary material). However, as all subsequent properties will be calculated from the same zone-center frequencies, we expect this to be intrinsically consistent with the remainder of the model.

It is worth highlighting that visualization of the cumulative total vibrational heat capacity, Fig. 3, with respect to increasing wavenumber provides another convenient way to visually locate Ω_{max}. Our assumption for Ω_{max} requires *g*(*ω*) = 0, which is given by a plateau in the cumulative heat capacity curves, Fig. 3(b). If the first such plateau is taken to be indicative of Ω_{max}—denoting the first break in *g*(*ω*)—reasonable agreement with our MC approach is obtained, Table II. In many cases, however, a small plateau is observed at wavenumbers well below predicted values of Ω_{max} by the MC method. To understand the sensitivity of our up-pumping model to the choice of Ω_{max}, we will consider the potential of these lower bounds; note, however, that the lower values generally do not bound the number of modes, which are amalgamated by the MC, EV-MC, or CoM approaches (see the supplementary material, Table II).

### Predicting impact sensitivity

A critical aspect of the up-pumping model, described above, pertains to the normalization of the phonon density of states, *g*(*ω*). Although a count of *g*(*ω*) is provided by the DFT simulations, it accounts only for pathways that reside at **k** = 0 or from an arbitrary number of **k** points. To reflect better the true relative number of vibrational pathways present in each system, we apply a normalization of $\u222b\Omega maxd\omega g\omega =Z6+YV0$. A solid with volume *V* and containing *N* unit cells comprises *N* wave vectors. Hence, in a region **k** + d**k**, the number of wave vectors [described by Eq. (8)] differs between systems only by the volume of the unit cell,

It, therefore, follows that normalizing by the unit cell volume offers a direct route for comparing the relative number of transfer pathways available across different EM crystalline solids.

To reflect the relative amount of energy transferred to each EM upon impact, we assume that the compressibility factor and Grüneisen parameters expressed in Eq. (1) are constant across all EMs (see the supplementary material). This is generally a good approximation based on experimentally available data.^{60} For the purpose of this work, we take as a model a value of *V*∖*V*_{0} ≈ 0.8, with a value for *γ* ≈ 4. This results in a final temperature *T*_{f} = 953.6 K and hence, a total bulk heating Δ*T* = 655.6 K. For each material, the corresponding initial quasi-temperature of the phonon path, *ϕ*(0)_{ph}, is obtained by Eq. (7) and is, hence, directly proportional to *C*_{tot}/*C*_{ph}, Table III (and the supplementary material).

. | . | EV-MC/CoM . | HC Model . | ||||
---|---|---|---|---|---|---|---|

. | C_{tot}
. | C_{ph}
. | C_{tot}/C_{ph}
. | ϕ(0)_{ph} (K)
. | C_{ph}
. | C_{tot}/C_{ph}
. | ϕ(0)_{ph} (K)
. |

HNB | 585.85 | 112.20 | 5.22 | 3423.3 | 112.20 | 5.22 | 3423.3 |

ϵ-CL-20 | 892.89 | 178.69 | 5.00 | 3278.0 | 157.92 | 5.65 | 3704.1 |

β-HMX | 685.48 | 128.82 | 5.32 | 3487.8 | 62.33 | 11.00 | 7211.6 |

o-TNT | 521.14 | 122.59 | 4.25 | 2786.3 | 89.35 | 5.83 | 3822.1 |

α-FOX-7 | 342.76 | 68.56 | 5.00 | 3278.0 | 68.56 | 5.00 | 3278.0 |

α-NTO | 255.45 | 63.33 | 4.03 | 2642.1 | 55.03 | 4.64 | 3042.0 |

TATB | 585.78 | 87.26 | 6.71 | 4399.1 | 87.26 | 6.71 | 4399.1 |

. | . | EV-MC/CoM . | HC Model . | ||||
---|---|---|---|---|---|---|---|

. | C_{tot}
. | C_{ph}
. | C_{tot}/C_{ph}
. | ϕ(0)_{ph} (K)
. | C_{ph}
. | C_{tot}/C_{ph}
. | ϕ(0)_{ph} (K)
. |

HNB | 585.85 | 112.20 | 5.22 | 3423.3 | 112.20 | 5.22 | 3423.3 |

ϵ-CL-20 | 892.89 | 178.69 | 5.00 | 3278.0 | 157.92 | 5.65 | 3704.1 |

β-HMX | 685.48 | 128.82 | 5.32 | 3487.8 | 62.33 | 11.00 | 7211.6 |

o-TNT | 521.14 | 122.59 | 4.25 | 2786.3 | 89.35 | 5.83 | 3822.1 |

α-FOX-7 | 342.76 | 68.56 | 5.00 | 3278.0 | 68.56 | 5.00 | 3278.0 |

α-NTO | 255.45 | 63.33 | 4.03 | 2642.1 | 55.03 | 4.64 | 3042.0 |

TATB | 585.78 | 87.26 | 6.71 | 4399.1 | 87.26 | 6.71 | 4399.1 |

Following normalization of the phonon density of states, our up-pumping model follows the two-step process described in Eq. (6). At each step of the up-pumping calculation, the up-pumped states are dispersed over the *Z* molecules, reflecting the localization of the vibrational energy with increasing frequency. A shock temperature (Table III) is initially applied to the underlying phonon DOS (pre-set with populations corresponding to 300 K), yielding a visibly “hot” phonon bath, Fig. 4(a). Following from Eq. (6a), this excess energy is propagated onto the states residing at Ω_{max} < *ω* < 2Ω_{max}, Fig. 4(b). Finally, Eq. (6b) allows for the final up-pumping of density to a maximum of *ω* < 3Ω_{max}, Fig. 4(c).

Propagation of vibrational density in this way provides a measure of the total excitation of the system. However, to yield initiation, the vibrational energy must induce an electronic response. This process is known as dynamical metallization, wherein heightened vibrational excitation of the target frequency induces structural changes that result in closing of the electronic band gap.^{37} In simple systems, such as metal azides, there exists a well-defined “target” vibrational mode that clearly facilitates the necessary metallization phenomenon. When considering more complex molecules, as done here, no single vibrational mode can be readily identified as being responsible for initiation. Instead, we adopt the arguments of RRKM kinetic theory of unimolecular decomposition. Noting that for large molecules the two-phonon density of states is approximately continuous over the range Ω_{max} − 3Ω_{max} (see the supplementary material), all up-pumped energy into this region can rapidly redistribute into the underlying fundamental modes present in the original density of states. Under this assumption, the probability of unimolecular decomposition is proportional to the total amount of up-pumped vibrational energy. Correspondingly, we take the integrated density of activated fundamental states as being indicative of the impact sensitivity, Fig. 5.

The results in Fig. 5(a) clearly indicate that our up-pumping model using the combined EV-MC or CoM-Ω_{max} values captures the relative experimental IS ordering for our EM test set well. Although negligible vibrational population is up-converted in the insensitive compounds *α*-NTO and TATB, the sensitive compounds HNB and *ϵ*-CL-20 exhibit large up-pumping contributions. Moreover, the ordering of the intermediate compounds *β*-HMX, o-TNT, and *α*-FOX-7 is also correctly predicted. The correct ordering of sensitive compounds HNB and *ϵ*-CL-20—which differ by only 0.25 Nm—is correctly captured, although a minor misordering of o-TNT and *α*-FOX-7 (ΔIS *ca* 1 Nm) is observed. We note, however, that such small differences in IS values are not generally meaningful and are well within the errors of the experimental testing approach. Even the use of the HC-Ω_{max} values provides a reasonable ordering outcome [Fig. 5(b)], with the most notable change being a reduced resolution in the relative predicted sensitivities for *β*-HMX and *ϵ*-CL-20. Despite these minor losses in the predictive capacity of the model, the overall classification of the test materials is, however, retained.

Our data, therefore, suggest that while the choice of Ω_{max} is undoubtedly important in predicting IS, it should be thought of as an indirect parameter only. From a materials design perspective, this also suggests that the initial material excitation (i.e., the value of *ϕ*_{ph}) is not directly indicative of the impact reactivity of the EMs tested here. This is further evident from the tabulated values of *ϕ*_{ph} in Table III and a lack of direct correlation with the impact sensitivity response. Correspondingly, features that directly influence the bulk and phonon heat capacities [*cf.* Eq. (3)] are unlikely to provide sensitive design targets for new EMs. For example, a low phonon shock quasi-temperature is expected for systems comprising simultaneously a low total heat capacity (induced by a small number of molecular degrees of freedom) and a high phonon heat capacity (induced by a high *Z* count in the primitive unit cell and a large number of *Y* amalgamated molecular vibrational modes). While these variables are undoubtedly important in understanding the vibrational up-pumping process, we must look into other areas to influence the structure/property relationship.

The predictions obtained from our model are instead dominated by the underlying vibrational structure. In particular, the model is dominated by the availability of vibrational states in the range Ω_{max} < *ω* < 2Ω_{max},^{35} i.e., the region that captures directly the energy up-pumped from the “hot” phonon bath, Eq. (6a). This dominant aspect of the up-pumping model has been noted previously in the literature.^{38} The up-pumped density in the region 2Ω_{max} < *ω* < 3Ω_{max} captures the second stage of the up-converted phonon energy, Eq. (6b), and also plays a crucial role in determining the predicted sensitivity ordering. It is, therefore, likely that influencing the distribution of vibrational states within these windows will provide the design-space for new energetic materials with tailored impact sensitivity responses.

The range of vibrational modes driving the phonon up-pumping model spans from ∼200 ± 50 cm^{−1} to 600 ± 150 cm^{−1}. This is the vibrational window determined largely by the flexibility of the energetic molecule, including angle bends, ring deformations, and other similar molecular distortions. Correspondingly, it is reasonable to expect that molecular flexibility will be indicative of the vibrational density in this critical window and hence, of the overall sensitivity of the material. To this end, we calculated two established flexibility indices: (1) the rotatable bond count (i.e., the number of single bonds that are not in a ring and that are bound to a non-terminal heavy atom) and (2) the Kier Molecular Flexibility (KMF) index given by^{61}

where $K\alpha 1$ and $K\alpha 2$ represent the first- and second-order shape indices [i.e., molecular graph parameters indicative of the atom (vertex) count, vertex cycles, and the degree of branching] and *N* is the number of vertex points in the molecule. Both molecular flexibility indices yield very similar data to each other and, with a few anomalous data points, provide a generally good indication of the relative sensitivity of the EMs in our test set, Fig. 6. The flexibility indices of TATB are on par with TNT, despite a pronounced difference in their IS response. Similarly, CL-20 exhibits a KMF that is nearly two orders of magnitude above the remaining EMs, although overall, it is correctly predicted to be highly sensitive. It follows that, designing energetic molecules with a high degree of molecular flexibility is likely to be a strong indication of overall sensitivity. This is consistent with our up-pumping model, where this flexibility is encoded into the molecular vibrations that fall into the Ω_{max} < *ω* < 3Ω_{max} window. This same flexibility can be influenced by the strength of interatomic interactions, which can hinder rotations and modify charge densities and bond dissociation energies. Thus, tuning intermolecular interactions through polymorphism or multi-component crystallization can also be expected to tune reactivity to mechanical impact.

It is worth considering briefly the effect of using only the zone-center phonon density of states to compute vibrational up-pumping sensitivity indicators. For comparison, the up-pumping IS indicators based on the full phonon dispersion curves for a subset of the EMs reported here are included in the supplementary material. The full phonon dispersion curves generally possess broader and more highly occupied phonon baths. Negligible differences between the datasets are expected beyond Ω_{max} as the mode behavior switches from lattice to molecular. Importantly, owing to minimal dispersion in the top phonon bath vibrational mode, values of Ω_{max} are not markedly affected. Numerical differences in our IS sensitivity indicators—largely due to changes in the phonon bath region—are, therefore, expected between the datasets and are in fact observed. Moreover, the difference between these data sets is system-dependent, reflecting the relative magnitude of phonon dispersion across our test systems. Critically, the zone-center and full phonon dispersion curve datasets produce the same unambiguous impact sensitivity ranking, meaning that fair ranking of predicted impact sensitivities appears sensible so long as the simulated phonon density of states have been obtained in a consistent way. As such, our work demonstrates how the use of the zone-center frequencies, which are substantially quicker to generate than the full dispersion curves, should be considered sufficient for a computational screening study.

## SUMMARY OF UP-PUMPING APPROACH

In summary, the proposed up-pumping model for IS prediction follows five principle steps, which are expected to be easily automated: (1) Following simulation of zone center vibrational frequencies, we elect to impose a Gaussian broadening of 5 cm^{−1} to approximate phonon dispersion. The Gaussian broadened spectrum is normalized by the total number of wave vectors in the material via Eq. (8). This provides the underlying vibrational spectrum upon which up-pumping occurs. (2) The initial conditions for up-pumping require the definition of a phonon bath, with upper bound Ω_{max}. Our results suggest that Ω_{max} can be readily obtained by analyzing the displacement of the molecular center of mass (CoM) for each simulated eigenvector. This offers a clear indication of the frequency that bounds all external phonon modes. (3) The quasi-temperature of the phonon bath is subsequently obtained by considering adiabatic compression, Eq. (1), in which all EMs are here assumed to have identical Grüneisen coefficients and compressibility. Correspondingly, the final bulk temperature of each material is equivalent. By calculating the total and phonon heat capacities, Eq. (7), the material-specific phonon quasi-temperature is obtained via *C*_{tot}/*C*_{ph}. (4) This excess phonon population is up-pumped via Eqs. (6a) and (6b). (5) Finally, as our approach is limited to first order anharmonic effects, we take the total quantity of the up-pumped population within the region Ω_{max} < *ω* < 3Ω_{max} as being indicative of the relative IS of the material.

## CONCLUSIONS

The search for new EMs that offer the correct balance of power, safety, and compliance with environmental legislation is an inherently hazardous endeavor. Herein, we have demonstrated that impact sensitivity, one of the key safety metrics for EMs, is amenable to prediction via a two-stage vibrational up-pumping model based on the Brillouin zone-center vibrational density of states. Our model correctly predicts the relative impact sensitivities for a test set of molecular crystals that are both structurally and energetically diverse and are some of the most well-known EMs in commercial and academic research fields. Moreover, this is achieved without resort to any fitting parameters, meaning that it can be applied to any energetic material, provided the crystal structure is known. As a critical decision in implementing the up-pumping approach, we have explored a series of possible selections for Ω_{max}, namely, a mode counting approach (MC), its combination with manual eigenvector analysis (EV-MC), a related approach based on screening displacements of the molecular center of mass (CoM), and analysis of cumulative heat capacities (HC). Although we expect the EV-MC approach to be most physically meaningful, it requires significant manual effort and is, therefore, not amenable to automated screening. The CoM approach provides the same values of Ω_{max} as EV-MC and is readily implemented for automated and high-throughput screening. We, therefore, expect the CoM methodology to become useful for screening libraries of potential EMs.

The up-pumping model goes beyond simply reproducing known experimental data. Having a predictive tool allows exploration of structure/property relationships at the most fundamental of levels and thereby provides a feedback loop to rationalize the material design with tailored energetic properties. Here, the distribution of molecular modes that reside within the 1–3Ω_{max} range hold the key to designing a material that is classed as either a primary (highly impact sensitive) or a secondary (lower impact sensitive) energetic. It, therefore, follows that molecular design, through the introduction (or removal) of low-energy amalgamated vibrational modes, and studies on polymorphism and multi-component crystallization, to increase (or reduce) the strengths of intermolecular interactions and charge distributions through the lattice, all have a part to play in shifting the pattern and distribution of low energy molecular vibrations. This, in turn, affects the efficiency of the molecular crystal to adsorb the energy from the mechanical impact. The vibrational up-pumping model draws a direct line from the IS prediction models based on the electronic structure through a physical model of energy transfer and bond dissociation. It is clear that having reliable computational modeling tools, to act in a screening capacity in the search for new energetic materials, is an important step forward in the field of energetic material research.

## SUPPLEMENTARY MATERIAL

See the supplementary material for the simulated crystallographic structures, computed vibrational spectra, and further details on the calculation of solid-state heat capacities and vibrational up-pumping processes.

## ACKNOWLEDGMENTS

The authors are grateful to the Edinburgh Computing and Data Facilities (ECDF), BAM IT, and the UK Materials and Molecular Modeling Hub, which was partially funded by EPSRC (Grant No. EP/P020194/1), for computing resources. This work has been performed under Project No. HPC-EUROPA3 (No. INFRAIA-2016-1-730897), with the support of the EC Research Innovation Action under the H2020 Program; in particular, A.A.L.M. gratefully acknowledges technical support provided by the EPCC, Edinburgh, UK. This material is based on the work supported by the Air Force Office of Scientific Research under Award No. FA8655-20-1-7000.

## DATA AVAILABILITY

The data that support the findings of this study are available within the article and its supplementary material.