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.

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 azides26 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

FIG. 1.

Schematic representation of molecular structures of the EMs used in this study.

FIG. 1.

Schematic representation of molecular structures of the EMs used in this study.

Close modal

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 excitation36 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 Bernstein38 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.

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/V0 and the Grüneisen parameter, γ. The final temperature Tf that is reached by an adiabatically compressed solid at initial temperature T0 with initial volume V0 can be described according to40 

TfT0=VV0γ.
(1)

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

H=T0TfdTC(T).
(2)

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 Cph such that

CphϕphT0=T0TfdTC(T),
(3)

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,

H=Hi+HL+Hc,
(4)

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

HcΣi,j,kVi,j,k(3)δΩiωjωk,

where

Vi,j,k(3)=3Euiujuk.
(5)

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 evidence42 and semi-empirical calculations43 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 developments33 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,

τ1=ρ2(ω=ωT/2)×[nphnd],
(6a)
τ2=ρ(2)×[nph+ndnT].
(6b)

Here, it is assumed that τ1 transfers all excess population energy (denoted by the Bose–Einstein populations, n) from the phonon states (nph) into the doorway region (nd), thereby re-establishing equilibrium populations of nph prior to further up-pumping. The “hot” nd states subsequently up-pump to the target modes (nT) 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 (nT) 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.

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 

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 grid47 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.

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 size53 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.

TABLE I.

Sensitivity and selected crystallographic parameters for the EM test set. References are provided for the reported impact sensitivity values.

EMIS (Nm)Vexp3)Vcalc3)ΔV (%)Space groupReferences
HNB 2.75 581.49 585.25 +0.6 I2/c 55  
ϵ-CL-20 3–4.25 1397.25 1461.01 +4.6 P21/c 56  
β-HMX 6.5–8 501.37 516.01 +2.9 P21/c 56 and 57  
o-TNT 24 1823.55 1839.64 +0.8 Pca21 58  
α-FOX-7 25.2–31.8 522.33 514.32 −1.5 P21/n 58 and 59  
α-NTO 72.75 902.06 927.22 +2.8 P1¯ 57  
TATB 122.5 425.25 434.69 +2.2 P1¯ 57  
EMIS (Nm)Vexp3)Vcalc3)ΔV (%)Space groupReferences
HNB 2.75 581.49 585.25 +0.6 I2/c 55  
ϵ-CL-20 3–4.25 1397.25 1461.01 +4.6 P21/c 56  
β-HMX 6.5–8 501.37 516.01 +2.9 P21/c 56 and 57  
o-TNT 24 1823.55 1839.64 +0.8 Pca21 58  
α-FOX-7 25.2–31.8 522.33 514.32 −1.5 P21/n 58 and 59  
α-NTO 72.75 902.06 927.22 +2.8 P1¯ 57  
TATB 122.5 425.25 434.69 +2.2 P1¯ 57  
TABLE II.

Estimation of Ωmax based on the mode-counting (MC), eigenvector mode-counting (EV-MC), center of mass (CoM), and heat capacity (HC) approaches. The number of molecules in the unit cell, Z, number of atoms per molecule, N, and the number of amalgamated modes per molecule, Y, are given.

Ωmax/cm−1
EMZ/unit cellN/moleculeYMCEV-MCCoMHCa
HNB 24 215 140 140 145 
ϵ-CL-20 36 16 217 217 217 180/222 
β-HMX 28 11 237 191 191 106/196 
o-TNT 21 212 212 212 160/217 
α-FOX-7 14 169 169 169 174 
α-NTO 11 251 199 199 166/204 
TATB 24 140 140 140 145 
Ωmax/cm−1
EMZ/unit cellN/moleculeYMCEV-MCCoMHCa
HNB 24 215 140 140 145 
ϵ-CL-20 36 16 217 217 217 180/222 
β-HMX 28 11 237 191 191 106/196 
o-TNT 21 212 212 212 160/217 
α-FOX-7 14 169 169 169 174 
α-NTO 11 251 199 199 166/204 
TATB 24 140 140 140 145 
a

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

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., NO2 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.

FIG. 2.

Phonon density of states for EMs in the test set. Ωmax as identified by the EV-MC and CoM methods is indicated with a dashed vertical line, and the gas phase frequencies below this point are shown as blue vertical lines.

FIG. 2.

Phonon density of states for EMs in the test set. Ωmax as identified by the EV-MC and CoM methods is indicated with a dashed vertical line, and the gas phase frequencies below this point are shown as blue vertical lines.

Close modal

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, Nint. Owing to the symmetry splitting of molecular modes in the solid state, this yields Y = ZNint 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 3N 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 Ctot/Cph. The high temperature limit of the bulk and phonon heat capacities can then be determined according to

cv=0ΩmaxnTωgωdω,
(7)

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).

FIG. 3.

Calculation of the vibrational heat capacity for TATB based on the phonon density of states, according to Eq. (7). (a) Overlay of the cumulative vibrational heat capacity and phonon density of states. (b) Identification of Ωmax based on the plateau in the cumulative heat capacity.

FIG. 3.

Calculation of the vibrational heat capacity for TATB based on the phonon density of states, according to Eq. (7). (a) Overlay of the cumulative vibrational heat capacity and phonon density of states. (b) Identification of Ωmax based on the plateau in the cumulative heat capacity.

Close modal

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 Ωmaxdωgω=Z6+YV0. A solid with volume V and containing N unit cells comprises N wave vectors. Hence, in a region k + dk, the number of wave vectors [described by Eq. (8)] differs between systems only by the volume of the unit cell,

gkdk=V2π34πk2dk.
(8)

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 VV0 ≈ 0.8, with a value for γ ≈ 4. This results in a final temperature Tf = 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 Ctot/Cph, Table III (and the supplementary material).

TABLE III.

Predicted total and phonon-bath only heat capacities and phonon bath shock temperature heating used in the IS prediction model, calculated based on EV-MC and HC derived values of Ωmax (shown in Table II; lowest values employed in the HC-model to explore the sensitivity of up-pumping to Ωmax).

EV-MC/CoMHC Model
CtotCphCtot/Cphϕ(0)ph (K)CphCtot/Cphϕ(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/CoMHC Model
CtotCphCtot/Cphϕ(0)ph (K)CphCtot/Cphϕ(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).

FIG. 4.

Key steps in the two-step up-pumping model for TATB. (a) The hot phonon bath and thermally equilibrated internal modes, (b) hot doorway states obtained by up-pumping through Eq. (6a), and (C) up-pumped density obtained via Eq. (6b) and projected onto fundamental bands. In each case, the phonon density of states g(ω) is shown in black, with the up-pumping effects on the populations shown in blue. In (b) and (c), the vertical dotted lines denote Ωmax and 3Ωmax.

FIG. 4.

Key steps in the two-step up-pumping model for TATB. (a) The hot phonon bath and thermally equilibrated internal modes, (b) hot doorway states obtained by up-pumping through Eq. (6a), and (C) up-pumped density obtained via Eq. (6b) and projected onto fundamental bands. In each case, the phonon density of states g(ω) is shown in black, with the up-pumping effects on the populations shown in blue. In (b) and (c), the vertical dotted lines denote Ωmax and 3Ωmax.

Close modal

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.

FIG. 5.

Predicted relative sensitivity using the vibrational up-pumping model. (a) Based on EV-MC/CoM-Ωmax values and (b) based on lowest HC-Ωmax values. Further details are available in the supplementary material, along with predicted sensitivities based on the MC-Ωmax values.

FIG. 5.

Predicted relative sensitivity using the vibrational up-pumping model. (a) Based on EV-MC/CoM-Ωmax values and (b) based on lowest HC-Ωmax values. Further details are available in the supplementary material, along with predicted sensitivities based on the MC-Ωmax values.

Close modal

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 by61 

KMF=Kα1Kα2N,
(9)

where Kα1 and Kα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.

FIG. 6.

Comparison of molecular flexibility parameters against impact sensitivity ordering, listed (left to right) in the order of decreasing impact sensitivity.

FIG. 6.

Comparison of molecular flexibility parameters against impact sensitivity ordering, listed (left to right) in the order of decreasing impact sensitivity.

Close modal

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.

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 Ctot/Cph. (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.

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.

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.

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.

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

1.
P. A.
Davies
,
J. Hazard. Mater.
38
,
75
(
1994
).
2.
Chemical Rocket Propulsion: A Comprehensive Survey of Energetic Materials
, edited by
A. S.
Cumming
(
Springer International Publishing
,
Cham
,
2017
), pp.
727
742
.
3.
United Nations, Committee of Experts on the Transport of Dangerous Goods, United Nations, and Committee of Experts on the Transport of Dangerous Goods and on the Globally Harmonized System of Classification and Labelling of Chemicals, Recommendations on the Transport of Dangerous Goods: Model Regulations,
2019
.
4.
Recommendations on the Transport of Dangerous Goods—Manual of Tests and Criteria Hauptbd
, revised edition, edited by
V.
Nationen
(
United Nations
,
New York
,
2009
), Vol. 5.
5.
R.
Tsyshevsky
,
O.
Sharia
, and
M.
Kuklja
,
Molecules
21
,
236
(
2016
).
6.
P.
Politzer
and
J. S.
Murray
,
Advances in Quantum Chemistry
(
Elsevier
,
2014
), pp.
1
30
.
7.
B. M.
Rice
,
S.
Sahu
, and
F. J.
Owens
,
J. Mol. Struct.: THEOCHEM
583
,
69
(
2002
).
8.
T. L.
Jensen
,
J. F.
Moxnes
,
E.
Unneberg
, and
D.
Christensen
,
J. Mol. Model.
26
,
65
(
2020
).
9.
D.
Mathieu
,
Propellants, Explos., Pyrotech.
45
,
966
973
(
2020
).
10.
M. J.
Kamlet
and
H. G.
Adolph
,
Propellants, Explos., Pyrotech.
4
,
30
(
1979
).
11.
V. W.
Manner
,
D. N.
Preston
,
C. J.
Snyder
,
D. M.
Dattelbaum
, and
B. C.
Tappan
,
AIP Conf. Proc.
1793
,
040036
(
2017
).
12.
M.
Jungová
,
S.
Zeman
, and
Q.-L.
Yan
,
Cent. Eur. J. Energ. Mater.
11
,
383
393
(
2014
), ISSN: 2353–1843.
13.
Z.
Jun
,
C.
Xin-lu
,
H.
Bi
, and
Y.
Xiang-dong
,
Struct. Chem.
17
,
501
(
2006
).
14.
J.
Mullay
,
Propellants, Explos., Pyrotech.
12
,
60
(
1987
).
15.
J. S.
Murray
,
M. C.
Concha
, and
P.
Politzer
,
Mol. Phys.
107
,
89
(
2009
).
16.
W. J.
Moodie
,
Safety Considerations in Grinding HMX
(
Mound Facility, Maimisburg, OH
,
1978
).
17.
P.
Soni
,
C.
Sarkar
,
R.
Tewari
, and
T. D.
Sharma
,
J. Energ. Mater.
29
,
261
(
2011
).
18.
R. K.
Weese
and
A. K.
Burnham
,
Propellants, Explos., Pyrotech.
30
,
344
(
2005
).
19.
S. R.
Kennedy
and
C. R.
Pulham
, in
Monographs in Supramolecular Chemistry
, edited by
C. B.
Aakeröy
and
A. S.
Sinha
(
Royal Society of Chemistry
,
Cambridge
,
2018
), pp.
231
266
.
20.
C. B.
Aakeröy
,
T. K.
Wijethunga
, and
J.
Desper
,
Chem.-Eur. J.
21
,
11029
(
2015
).
21.
M.
Cartwright
and
J.
Wilkinson
,
Propellants, Explos., Pyrotech.
36
,
530
(
2011
).
22.
W.
Zhu
and
H.
Xiao
,
Struct. Chem.
21
,
657
(
2010
).
23.
S.
Zeman
,
N.
Liu
,
M.
Jungová
,
A. K.
Hussein
, and
Q.-L.
Yan
,
Def. Technol.
14
,
93
(
2018
).
24.
P.
Politzer
and
J. S.
Murray
,
Struct. Chem.
27
,
401
(
2016
).
25.
S. V.
Bondarchuk
,
CrystEngComm
20
,
5718
(
2018
).
26.
S. V.
Bondarchuk
,
New J. Chem.
43
,
1459
(
2019
).
27.
S. V.
Bondarchuk
,
J. Phys. Chem. A
122
,
5455
(
2018
).
28.
Y.
Ma
,
A.
Zhang
,
C.
Zhang
,
D.
Jiang
,
Y.
Zhu
, and
C.
Zhang
,
Cryst. Growth Des.
14
,
4703
(
2014
).
29.
J.
Zhang
,
L. A.
Mitchell
,
D. A.
Parrish
, and
J. M.
Shreeve
,
J. Am. Chem. Soc.
137
,
10532
(
2015
).
30.
C.
Zhang
,
X.
Wang
, and
H.
Huang
,
J. Am. Chem. Soc.
130
,
8359
(
2008
).
31.
R.
Bu
,
Y.
Xiong
, and
C.
Zhang
,
Cryst. Growth Des.
20
,
2824
(
2020
).
32.
Y.
Ma
,
A.
Zhang
,
X.
Xue
,
D.
Jiang
,
Y.
Zhu
, and
C.
Zhang
,
Cryst. Growth Des.
14
,
6101
(
2014
).
33.
D. D.
Dlott
and
M. D.
Fayer
,
J. Chem. Phys.
92
,
3798
(
1990
).
34.
A. A. L.
Michalchuk
,
P. T.
Fincham
,
P.
Portius
,
C. R.
Pulham
, and
C. A.
Morrison
,
J. Phys. Chem. C
122
,
19395
(
2018
).
35.
A. A. L.
Michalchuk
,
M.
Trestman
,
S.
Rudić
,
P.
Portius
,
P. T.
Fincham
,
C. R.
Pulham
, and
C. A.
Morrison
,
J. Mater. Chem. A
7
,
19539
(
2019
).
36.
A.
Tokmakoff
,
M. D.
Fayer
, and
D. D.
Dlott
,
J. Phys. Chem.
97
,
1901
(
1993
).
37.
A. A. L.
Michalchuk
,
S.
Rudić
,
C. R.
Pulham
, and
C. A.
Morrison
,
Phys. Chem. Chem. Phys.
20
,
29061
(
2018
).
38.
J.
Bernstein
,
J. Chem. Phys.
148
,
084502
(
2018
).
40.
A. I.
Kitaigorodsky
,
Molecular Crystals and Molecules
(
Academic Press
,
1973
).
41.
C. S.
Coffey
and
E. T.
Toton
,
J. Chem. Phys.
76
,
949
(
1982
).
42.
S. D.
McGrane
,
J.
Barber
, and
J.
Quenneville
,
J. Phys. Chem. A
109
,
9919
(
2005
).
43.
Ye
and
M.
Koshi
,
J. Phys. Chem. B
110
,
18515
(
2006
).
44.
F.
Neese
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
8
,
e1327
(
2018
).
45.
M. K.
Kesharwani
,
B.
Brauer
, and
J. M. L.
Martin
,
J. Phys. Chem. A
119
,
1701
(
2015
).
46.
S. J.
Clark
,
M. D.
Segall
,
C. J.
Pickard
,
P. J.
Hasnip
,
M. I. J.
Probert
,
K.
Refson
, and
M. C.
Payne
,
Z. Kristallogr. - Cryst. Mater.
220
,
567
(
2005
).
47.
H. J.
Monkhorst
and
J. D.
Pack
,
Phys. Rev. B
13
,
5188
(
1976
).
48.
S.
Hunter
,
P. L.
Coster
,
A. J.
Davidson
,
D. I. A.
Millar
,
S. F.
Parker
,
W. G.
Marshall
,
R. I.
Smith
,
C. A.
Morrison
, and
C. R.
Pulham
,
J. Phys. Chem. C
119
,
2322
(
2015
).
49.
S.
Hunter
,
T.
Sutinen
,
S. F.
Parker
,
C. A.
Morrison
,
D. M.
Williamson
,
S.
Thompson
,
P. J.
Gould
, and
C. R.
Pulham
,
J. Phys. Chem. C
117
,
8062
(
2013
).
50.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
51.
S.
Grimme
,
J. Comput. Chem.
27
,
1787
(
2006
).
52.
K.
Refson
,
P. R.
Tulip
, and
S. J.
Clark
,
Phys. Rev. B
73
,
155114
(
2006
).
53.
W. A.
Trzciński
,
S.
Cudziło
,
Z.
Chyłek
, and
L.
Szymańczyk
,
J. Energ. Mater.
31
,
72
(
2013
).
54.
S.
Konar
,
A. A. L.
Michalchuk
,
N.
Sen
,
C. L.
Bull
,
C. A.
Morrison
, and
C. R.
Pulham
,
J. Phys. Chem. C
123
,
26095
(
2019
).
55.
W. S.
Wilson
,
D. E.
Bliss
,
S. L.
Christian
, and
D. J.
Knight
,
Explosive Properties of Polynitroaromatics
(
Defense Technical Information Center
,
Fort Belvoir, VA
,
1990
).
56.
R. L.
Simpson
,
P. A.
Urtiew
,
D. L.
Ornellas
,
G. L.
Moody
,
K. J.
Scribner
, and
D. M.
Hoffman
,
Propellants, Explos., Pyrotech.
22
,
249
(
1997
).
57.
C. B.
Storm
,
J. R.
Stine
, and
J. F.
Kramer
, in
Chemistry and Physics of Energetic Materials
, edited by
S. N.
Bulusu
(
Springer Netherlands
,
Dordrecht
,
1990
), pp.
605
639
.
58.
B. M.
Rice
and
J. J.
Hare
,
J. Phys. Chem. A
106
,
1770
(
2002
).
59.
I. J.
Lochert
,
FOX-7—A New Insensitive Explosive
(
DSTO Aeronautical and Maritime Research Laboratory
,
Victoria, Australia
,
2001
).
60.
Static Compression of Energetic Materials
, edited by
S. M.
Peiris
and
G.
Piermarini
(
Springer
,
Berlin
,
2008
).
61.
L. B.
Kier
,
Quant. Struct.-Act. Relat.
8
,
221
(
1989
).

Supplementary Material