Towards Computational Screening for New Energetic Molecules: Calculation of Heat of Formation and Determination of Bond Strengths by Local Mode Analysis

Introduction

Energetic materials (explosives, propellants, pyrotechnics and gas generators, EMs) are characterised by their ability to rapidly convert chemical potential energy into kinetic energy. Safety is of paramount importance in this field and takes precedence over performance, such that only a limited number of EMs have found employment in civilian and military applications. Many of the well-established EMs, however, do not comply with increasing environmental and public health regulations (Cumming, 2017), which fuels the need to develop new non-toxic and environmentally benign EMs that do not compromise on safety or performance metrics.

Computational screening presents an attractive route in the new materials pipeline, as it offers a cost-effective way to assess candidates prior to synthesis. This is particularly desirable in the field of EMs where safety testing (such as impact, spark and friction sensitivity measurements) typically require gram-scale quantities of compounds to be synthesised. This is potentially very hazardous work when the material is novel and may initiate on the slightest perturbation.

Having access to reliable predictive models also opens up routes to the rational design of new EMs, by offering a path towards understanding structure-property relationships. Previous work in our group has focused on the ability to predict the impact sensitivities of EMs using first principles simulations, and our methods, which are based on a vibrational up-pumping model, have demonstrated success for a range of structurally diverse materials that exhibit widely varying sensitivities to initiation ( Michalchuk et al., 2018a; Michalchuk et al., 2019; Michalchuk et al., 2021). This predictive model highlighted the importance of low energy (ca. 200 ± 50 to 600 ± 150 cm⁻¹) molecular vibrational modes to channel (up-pump) the energy arising from the phonon scattering of the many low energy lattice mode vibrations which become vibrationally hot following a mechanically-induced shock event. Trapping this energy in the low energy molecular vibrations induces large amplitude vibrations that distort the molecular structure to the degree that electronic changes occur: band gaps narrow, electrons flow, and unstable species emerge all on the timescale of a molecular vibration (Michalchuk et al., 2018b). This marks the start of initiation. It therefore follows that crystal lattices with a high density of low-lying molecular vibrations will likely be sensitive to shock and impact-induced initiation. The reverse also holds true: crystal lattices that are vibrationally sparse in this region will likely be shock and impact-insensitive, and are thus safer to handle.

Predicting impact sensitivities is only one aspect of an EM computational screening process. Also of key importance is the stored energy in the molecule, which can be gleaned from the solid-state heat of formation energy, ∆H_f(s). This provides a route to calculate the heat of combustion, which in turn allows prediction of several parameters including the detonation pressure and the velocity and heat of explosion (Politzer et al., 2003), using thermo-chemical software methods such as CHEETAH (Fried and Souers, 1994) or EXPLO-5 (Suceska, 2004).

Multiple routes to calculating ∆H_f(s) have been previously proposed in the literature. Our purpose here is to signpost the current state of the art, and to show its application to EMs. The first step in calculating ∆H_f(s) begins with the gas-phase heat of formation, ∆H_f(g), and various methods have been proposed for this over the years. For instance, Benson’s group increment theory (BIGT) exploited experimental heats of formation for individual groups of atoms to develop group equivalence values for linear and branched alkanes and alkenes (Benson and Buss, 1958). While impressive at the time of inception, its application is limited to the types of molecules represented in the training set, and so it has limited scope beyond this area. More recently, quantitative structure-property relationship (QSPR) models have been realised as a powerful tool to explore the relationships that link molecular structure to material properties. Vatani et al. (2007) devised a new QSPR method for predicting heats of formation for over 1,000 organic molecules, covering almost all organic functional groups. However, transition-metal and main-group elements were not included in the training set, so whilst excellent results were obtained, widespread application is again limited.

Within the field of EMs isodesmic reaction equations (Hehre et al., 1970; Ponomarev and Takhistov, 1997) and the atom- and group-equivalence methods (Byrd and Rice, 2006) are commonly applied. The former relies on a reaction equation where the types of chemical bonds in the reactants and products are conserved, and the heats of formation are known for all other molecules in the reaction except the one unknown (Cramer, 2004). Any intrinsic error associated with calculating any particular chemical bond is thereby cancelled out, meaning that relatively low-level computational methods can be employed to give fairly accurate results (Hehre et al., 1970). This method has been employed in computational chemistry for over 50 years, with recent developments automating the generation of parts of the isodesmic reaction equation (Chan et al., 2020). However, as a general technique it has found little application beyond CHNO-containing molecules.

The atom-equivalence method developed by Byrd and Rice is an advance on BGIT and has been extensively applied to EMs (Singh et al., 2014a; Singh et al., 2014b; Axthammer et al., 2015; Piercey et al., 2016), but is again restricted to CHNO-containing molecules. Here atom equivalence energy values for the four atoms were determined through comparison of experimental heats of formation for molecules in a training set and their computationally derived molecular energies (Byrd and Rice, 2006). This method works well, reporting root mean square deviation of 12.6 kJ mol⁻¹ and is arguably more efficient than the isodesmic equation route, as it requires just the optimised energy of the molecule of interest expressed at a prescribed quantum mechanical model chemistry. A further advantage is the absence of any reliance on further experimental data. However, the continued application to CHNO-only molecules limits its application in a computational screening programme for new EMs which should have the flexibility to draw upon main-group and transition-metal elements.

Given the limitations of these two methods, herein we have pursued the use of PM7 to determine ∆H_f(g) (Stewart, 2013). This has been shown by Wan et al. to out-perform previous semi-empirical methods for a set of 142 organic molecules (Wan et al., 2020). Elioff et al. evaluated its capabilities compared to both the isodesmic and the group equivalence methods for nitrogen-containing organic molecules (Elioff et al., 2016). While the outcome showed that PM7 was the least accurate of the three (R² line of best fit against experimental data = 0.986, compared to 0.999 and 0.995 for the isodesmic and atom equivalence methods, respectively), it still performs relatively well, and has the advantage of being an easily deployed method that can be used for any molecule containing atoms from H−La and Lu−Bi. Fomin et al. (2017) tested the PM7 method alongside other semi-empirical methods, namely PM6, PM5, PM3, AM/1, and MNDO, for calculating ∆H_f(g) for copper and alkaline Earth metal complexes. They concluded that PM6 and PM7 both perform well, reporting R² values of 0.961 and 0.960, respectively. While these results are encouraging, further validation for the accurate prediction of ∆H_f(g) for a broader range of inorganic molecules would be welcome and will be provided here. The PM7 method also carries the advantage that no further calculations beyond a geometry optimisation are required, which renders it attractive as part of a high throughput study. Moreover, semi-empirical calculation methods have a wide application and user base, and are being continuously improved (Wan et al., 2020).

Converting ∆H_f(g) to ∆H_f(s) requires the addition of an intermolecular interaction energy term, as captured by the sublimation energy, ∆H_sub, or the lattice energy, ∆H_L (formally, ∆H_L = −∆H_sub – 2RT). For single-component molecular solids, Politzer et al. have developed a route to determine ∆H_sub from consideration of the electrostatic potential (ESP) (Politzer et al., 1997; Politzer et al., 1998). While they have applied this analysis to estimate many liquid, solid, and solution parameters dependent on non-covalent interactions (Politzer and Murray, 1998; Politzer and Murray, 2002) it is their relationship to calculate ∆H_sub which is most relevant here. The method was developed initially using a training set of 34 CHNO-containing molecules, and further parametrised by Byrd et al. using a training set of 38 CHNO-containing energetic molecules (Byrd and Rice, 2006). The ESP method tested here uses the parameters proposed by Byrd, as parameterisation was carried out at using the B3LYP/6-31G* computational model, rather than the HF/STO-5G* level originally used by Politzer. In addition, Byrd’s training set included functional groups more likely to be present in EMs (e.g., azides and nitro groups).

For salts, an attractive route to ∆H_f(s) from ∆H_f(g) is via calculation of ∆H_L using the method developed by Jenkins et al. (1999), Jenkins et al. (2002) that relies only on knowledge of the molecular (formula unit) volume and the stoichiometry of the salt. This method has been used by Gao et al. (2007) to compare the estimated values of ∆H_f(s) for 33 energetic salts to their respective experimental values, using the G2 method to calculate the ∆H_f(g) terms for the ions based on their proton or electron affinities. They reported R² = 0.983, with a maximum deviation of 158.5 kJ mol⁻¹.

Co-crystals are an important development in the field of EMs (Kennedy and Pulham, 2018). While acknowledging that predicting whether or not a stable co-crystal will form from its single-component species is in itself the greater challenge, having the ability to predict ∆H_f(s) for known materials is also important. This is particularly true for EM research, where creating materials that can store large amounts of chemical energy is an essential requirement. Previous attempts have met with limited success. For instance, Zhang et al. (2016), Bozkuş et al. (2019) both used an atomisation energy method (Curtiss et al., 1997) which can formally only be used to calculate the heat of formation of gas-phase species. Ma et al. (2017) used isodesmic equations to calculate ∆H_f(g) for the co-formers of a CL-20/MTNP co-crystal, and then predicted ∆H_sub using a relationship based on the melting point of the co-crystal. While this method is promising, its application is limited if the melting point is not known. Gavezzotti developed the PIXEL method to calculate ∆H_L as a sum of Coulombic, polarisation, dispersion and repulsion intermolecular energies in a crystal structure and has been shown to have the correlation R² = 0.845 between experimental and calculated ∆H_L for 154 organic crystal structures (Gavezzotti, 2011). The method was recently expanded to include parameterisation for transition metal complexes (Maloney et al., 2015). Another route to obtaining ∆H_L is through the more computationally demanding dispersion-corrected density functional theory (DFT-D) (Morrison and Siddick, 2003; Feng and Li, 2006; Brandenburg et al., 2015; Fan et al., 2017). For the specific example of co-crystals considered here, as it is uncommon to measure ∆H_f(s) we shall compute formation enthalpies using both PIXEL and DFT-D for direct comparison.

In addition to having reliable routes to predict ∆H_f(s), having knowledge of the strengths of the individual bonds within a molecule is also valuable information at the EM molecular design stage, as it can provide information on the first-stage initiation pathway. Historically, intramolecular bond strengths have been calculated through heterolytic bond cleavage reactions, but this proves problematic beyond the first bond breaking reaction, as each subsequent bond breakage reaction is performed on increasingly unstable molecular fragments (Cremer et al., 2000), or requires a separate bond breaking reaction for each bond to be investigated (Zou et al., 2012). This is particularly problematic for ring systems, where breaking one bond introduces additional strain in the remaining bonds, such that isolating one bond becomes an impossible task. Alternatively, computation of the bond force constants offers a direct route to determining the bond strengths of all bonds within a molecule without recourse to molecular fragmentation. However, the normal modes of molecular vibration, from which the force constants are extracted, are commonly a complex mix of bond stretching, angle bending and twisting motions, meaning that pure bond force constants can rarely be obtained. Recent developments by Konkoli and Cremer (Konkoli and Cremer, 1998a; Konkoli and Cremer, 1998b; Kraka et al., 2020) have allowed for the mass decoupling of the normal modes of vibration, to recast the eigenvectors onto a new set of modes, termed the local modes of vibration, that correlate directly with individual bond stretches, angle bends etc. Their work has shown that the resulting local-mode force constants thus obtained are a direct measure of bond strength (Kraka and Cremer, 2019). Thus performing local mode analysis across a broad range of molecules (both energetic and non-energetic) provides information on the relationship between bond length and bond strength of the most common chemical bonds in EMs, which has the potential to be utilised for molecular design.

Herein, a set of 20 CHNO-containing molecules and a further 31 inorganic molecules was constructed to benchmark the PM7 method against isodesmic equation reactions and the atom equivalence method to calculate ∆H_f(g). Additionally, methods for converting ∆H_f(g) to ∆H_f(s) using the methods proposed by Byrd, Jenkins, and by the PIXEL method/DFT-D for single-component molecular crystals, salts and co-crystals, respectively, were also pursued for 48 compounds. Local vibrational mode analysis has also been carried out on 30 molecules containing chemical bonds found in energetic molecules, to evaluate bond length/strength relationships and to ascertain the likely weakest bonds in energetic molecules. Finally, both parameters have been highlighted for their potential to be included in a computational screening program for new energetic materials.

Computational Methods

All optimization and vibrational frequency calculations were performed at the B3LYP/6-31G* level, as implemented in Gaussian16 (Frisch et al., 2016).

Gas-Phase Heat of Formation: Isodesmic Equations

Equations were devised to ensure the type of chemical bonds were conserved, and that the heats of formation of all other molecules in the equation were known (see Supplementary Table 2). The heat of reaction, ∆H_R, is then calculated according to Eq. 1.

$\Delta {\mathit{H}}_{\mathit{R}}=\Delta {\mathit{E}}_0+\Delta \mathit{ZPE}+\Delta {\mathit{H}}_{\mathit{T}}+\Delta \mathit{nRT}(1)$

Where ΔE₀ is the change in energy between the products and reactants, ΔZPE is the change in the zero-point energies between products and reactants, and ΔH_T is the thermal correction from 0 to 298 K. As the number of atoms remain constant in the reaction, ΔnRT equals zero. The calculated heat of reaction is then equated to Eq. 2. Assuming the molecule of interest is a reactant in the isodesmic equation then its heat of formation is calculated by subtracting the known heats of formation of the other reactants and products from ∆H_R.

$\Delta H_R={\displaystyle \sum \Delta H_{f\left(g\right)products}} +{\displaystyle \sum \Delta H_{f\left(g\right)reactants}} (2)$

Gas-Phase Heat of Formation: Atom Equivalence Method

Formation energies were determined according to Eq. 3, where E is the optimised energy of the molecule, n_j is the number of atoms of type j and e_j is the atom equivalence value of atom j, as determined by Byrd and Rice (Byrd and Rice, 2006).

$\Delta H_{f\left(g\right)}=E-{\displaystyle \sum n_j{e}_j} (3)$

Gas-Phase Heat of Formation: PM7

PM7 (Stewart, 2013) method was utilised as presented in Gaussian 16, using geometries optimised to global minima from the previously mentioned calculations for improved accuracy. For molecules containing third row and higher atoms, the SCF = YQC algorithm was used, as suggested in the PM7 documentation.

Solid-State Heat of Formation: Single-Component Solids

The solid heat of formation for single component materials were calculated using Eq. 4, where ∆H_f(g) was calculated using PM7 as described above, and ∆H_sub was calculated using the ESP method as described by Eq. 5.

$\Delta H_{f\left(s\right)}=\Delta H_{f\left(g\right)}-\text{Δ}{H}_{sub}(4)$

$\text{Δ}{H}_{sub}=a{\left(SA\right)}^2+b\sqrt{{\sigma }_{tot}^2\nu } +c(5)$

Where $a, b$ and $c$ are semi-empirically deduced fitting parameters proposed by Byrd et al. (Byrd and Rice, 2006) $SA$ is the surface area of the 0.001 electron.bohr⁻³ isosurface of the electrostatic potential of the molecule, ${\sigma }_{tot}^2$is the measure of variablity of electronic potential on the surface, and $\nu$ is the degree of balance between the positive and negative charges on the isosurface. The latter three parameters were calculated using Multiwfn (Lu and Chen, 2012).

Solid-State Heat of Formation: Salts

For the energetic salts ΔH_f(s) was calculated from Eq. 6, where ∆H_f(g) of the cations and anions were calculated using PM7 as described above.

$\Delta H_{f\left(s\right)}=\Delta H_{f\left(g\right)\left(cation\right)}+\Delta H_{f\left(g\right)\left(anion\right)}-\Delta H_L(6)$

Here, ∆H_L is expressed by Eq. 7, as proposed by Jenkins et al. (2002).

$\Delta H_L=U_{pot}+\left[p\left(n_m/2-2\right)+q\left(n_x/2-2\right)\right]RT(7)$

Where $n_m$ and $n_x$ are constants that depend on the nature of the ions, and are set to 3 for monoatomic ions, 5 for linear polyatomic ions and 6 for non-linear polyatomic ions. The variables $p$ and $q$ denote the relative charges of the respective ions. The term $U_{pot}$ denotes the lattice potential energy and in turn is defined by Eq. 8.

$U_{pot}=\alpha {\left(V_m\right)}^{\frac{1}{3}}+\beta (8)$

Where $V_m$denotes the molecular volume (V_cell/Z), in nm³, and the remaining coefficients $\alpha$ and$\beta$ are fitting terms provided by Jenkins et al. and which vary depending on the charge ratio of the salt.

Solid-State Heat of formation: Co-Crystals

For lattice enthalpies calculated using PIXEL (Gavezzotti, 2005), calculations were set up using MrPIXEL (Reeves et al., 2020) within the Mercury interface, distributed with the Cambridge Structural Database (CSD) (Macrae et al., 2008; Groom et al., 2016). Hydrogen atom positions were set to the CSD normalised positions. For DFT-D, geometry optimisation calculations were performed using CASTEP17 (Clark et al., 2005) using the Perdew-Burke-Ernzerhof functional (Perdew et al., 1996) with a plane-wave basis set with a cut-off energy of 900 eV, which demonstrated convergence to 1 meV.atom⁻¹. Norm-conserving pseudopotentials were used throughout, with a k-point spacing of 0.05 Å⁻¹. The Tkatchenko-Scheffler dispersion correction scheme was applied. Lattice energies were determined by comparing the optimized energy values for the crystal structure with the energy for individual co-formers, modelled as effectively gas phase by removing all but one of the co-former molecules from the optimized crystal structure and computing a single point energy value using the same computational model as applied to the co-crystal. In cases where a unit cell vector was short, such that interactions with the nearest neighbour replica may occur (taken to be <5 Å), the smallest unit cell vector was doubled to ensure zero interaction. This was the case for 64, 70 and 76. Lattice enthalpies were then calculated using Eq. 9.

$\Delta H_L=\frac{E_{cell}}{Z}-E_{co-former1}-E_{co-fomer2}(9)$

Where E_cell is the energy of the unit cell of the co-crystal and E_co-fomer1,2 is the energy of each of the co-formers modelled in the “gas phase” and Z is the number of molecular units in the co-crystal. The ∆H_f(g) terms, to convert the ∆H_L terms to ∆H_f(s), were calculated using PM7 as documented above.

Lattice enthalpies of individual co-formers were calculated by adding the thermodynamic correction shown in Eq. 10 to the sublimation enthalpies calculated using the ESP method described above.

$\Delta H_L= -\Delta H_{sub}-2RT(10)$

Local Mode Force Constants

These were calculated using LMODEA (Kraka et al., 2020), following geometry optimisation and vibrational frequency calculation, for the 31 CHNO-containing molecules listed in Supplementary Table 6.

Results and Discussion

Calculating Gas Phase Heats of Formation

To test the three different methods for calculating ∆H_f(g) molecules 1–20 (see Figure 1 and Supplementary Table 1), were considered. Molecules were chosen as they had reliably reported experimental ∆H_f(g) available in the literature, and their less complex nature meant that the construction of isodesmic reaction equations was relatively straightforward (reported in full in Supplementary Table 2). As the PM7 method permits inorganic molecules to be studied, a further set of molecules which comprised main group and transition metal elements (23–53) was also studied. Notably, examples include molecules containing lead, copper and halogens, which are commonly found in EMs.

FIGURE 1

FIGURE 1. Molecules, salts and co-crystals employed in heat of formation energy calculations.

While the majority of CHNO-containing molecules investigated here have already been reported by Elioff et al. (2016), our process differs in that all molecules were first optimised at the 6-31G*/B3LYP level, followed by a single-point energy evaluation at PM7. In addition, the data set includes a further 4 molecules, namely 1, 2, 4, and 19, which incorporates the important EMs PETN (1), TMETN (2) into the test set. Calculation of ∆H_f(g) for the 31 inorganic molecules (23–53) by PM7 has not been reported before.

When calculated values of ∆H_f(g) are compared against experimental data (Figure 2), the two largest outliers for the isodesmic equation method (Figure 2A and Supplementary Table 3) were PETN and TMETN (1 and 2 respectively), which deviate from the experimental values by 106.8 and 108.0 kJ mol⁻¹, respectively. Both data points were disregarded when carrying out linear fitting, which otherwise returned an R² value of 0.994, and the gradient of the fitted line, m = 0.989. The reasons for failure for data points 1 and 2 must rest with either the experimental formation energies and/or the geometry optimisations of PETN and TMETN, or the geometry optimisation and/or experimental formation energy of another molecule defined in their respective isodesmic equations. The atom equivalence method (shown in Figure 2B) performs considerably better for these two compounds, which suggests that the experimental formation energies and the geometry optimizations, for both PETN and TMETN, are reliable. The isodesmic reactions constructed for these two molecules both include C(CH₃)₄, which is absent in all other isodesmic equations constructed for the test set (see Supplementary Table 2). Careful checking of the simulated geometry, to ensure the configuration obtained refers to the global minimum, suggests that the error most likely rests with the experimental heat of formation for C(CH₃)₄. This highlights a fundamental weakness with the isodesmic equation route, in that any error with any one term in the equation will render the calculated ∆H_f(g) for the target molecule as unreliable.

FIGURE 2

FIGURE 2. Showing calculated vs. experimental gas phase heats of formation for the (A) isodesmic, (B) atom equivalence, (C) PM7 for 1–20 and (D) PM7 for 23–53. Data points omitted from the lines of best fit are shown in pink.

Using the atom equivalence method, no outliers where identified, showing the strength of this method for CHNO-containing molecules. This gave an R² value of 0.994 (m = 0.943), indicating that the atom equivalence method performs just as well as the isodesmic method.

Turning to the PM7 results (Figure 2C), PETN (1) again represents the largest deviation, with the semi-empirical method over-estimating the experimental value by 91.2 kJ mol⁻¹. RDX (17) also appears to be less reliably calculated by this method, compared to the other two routes. Discarding 1 and 17 from the test set gives a line of best fit with R² = 0.993 (m = 1.083), suggesting overall that this method provides a similar level of accuracy compared to the isodesmic reaction and atom equivalence methods. The correlation with experimental data here is better than that reported by Elioff et al. (2016) (R² = 0.986) and considering the high overlap of molecules in test sets used, the improvement is likely due to the B3LYP/6-31G* geometry optimization step, as Elioff et al. relied on the PM7 method for structure optimisation as well as energy calculation. For the inorganic molecule list, 23–53, which were tackled by PM7 only, the R² value increased slightly to 0.995, with the gradient of the line of best fit improved considerably [m = 0.983; omitting the largest outliers ZrF₄ (24) and TeF₆ (26), which deviate from experimental values by 186 and 203 kJ mol⁻¹, respectively], suggesting its performance is just as reliable for inorganic molecules as it is for organic molecules. It is unclear why ZrF₄ and TeF₆ deviate so much from the expected values. Other related compounds, 25, 27–29 and 34, have been calculated accurately (although 26 is the only Te-containing compound and Zr is represented just twice in the data set). The possibility of experimental inaccuracy also cannot be ruled out.

Calculating Solid Heats of Formation

To probe the conversion of ∆H_f(g) to ∆H_f(s), 23 of the single-component EMs were selected from the 1–53 test set for which experimentally determined ∆H_f(s) were available (see Figure 3A and Supplementary Table 1), and further pursued using the ESP method. This includes 14 CHNO-containing molecules and 9 inorganic compounds. Fourteen salts (54–67 in Figure 1) were similarly pursued using Jenkins’ method, while 11 co-crystals (68–78 in Figure 1) were explored using PIXEL and DFT-D.

FIGURE 3

FIGURE 3. Showing the calculated solid heats of formation for (A) single component crystals, (B) salts and (C) co-crystals. Data points omitted from the lines of best fit are shown in pink.

For the single-component crystals, ∆H_f(s) calculated using the ESP method offered a line of best fit of R² = 0.995, m = 1.100 [omitting one point that lay off the line, Mn(CO)₃Cp (35) by 220 kJ mol⁻¹, see Figure 3A], indicating a good predictive result has been obtained across a broad spectrum of compounds and broad range of energies.

For the EM salts (54–67), application of Jenkins’ method gives a line of best fit through the simulated values (see Figure 3B) with R² = 0.997 (m = 0.955), and a maximum deviation from experimental values of 80.5 kJ mol⁻¹ for FOX-12 (55). This is an improvement in the correlation reported by Gao et al. (2007) who reported an R² of 0.984 for a set of 33 inorganic and CHNO EM salts, of which 19 were inorganic and the remainder CHNO. While acknowledging that the improved result obtained here could be due to the smaller test set employed in this work, it could also indicate that the PM7 method, used to calculate the ∆H_f(g) terms for the constituent ions, offered an improvement over the approach adopted by Gao, who relied on isodesmic equations. This latter approach is likely to be particularly problematic here as experimental formation energies of ions are required. Our data set shares two data points with Gao’s (salts 56 and 61); closer inspection of these predictions shows that both ∆H_f(s) values were calculated more accurately here, with 56 deviating from the experimental values by 17 kJ mol⁻¹ (Gao’s prediction differed by 24.7 kJ mol⁻¹), whereas 61 deviated by only 0.9 kJ mol⁻¹ (Gao by 63.5 kJ mol⁻¹). This would therefore appear to support further the use of PM7 to calculate the ∆H_f(g) terms. There are two outliers in this fitting – lead azide (64) and DBX-1 (65), which deviate from their literature values by 222 and 145 kJ mol⁻¹, respectively. We note that both exist as extended coordination complexes in the solid state, which may contribute towards their poor prediction by Jenkins’ method, which was formally devised for salts.

It should be noted that TKX-50 (62) has two values for ∆H_f(s) reported in the literature. The most widely reported is 446.6 kJ mol⁻¹ (Sabatini and Oyler, 2015), which was the value calculated using Jenkins’ method, with the ∆H_f(g) values for the constituent ions calculated using the CBS-4M atomisation method. However, Sinditskii et al. (2015) have argued that this value is questionable when compared to the sum of the enthalpies of formation for the individual components of TKX-50 and when compared with typical heats of reaction between acids and bases. They performed bomb calorimetry experiments and determined ∆H_f(s) to be 111 ± 16 kJ mol⁻¹, far lower than the widely reported calculated value. Our computed value, also derived from the Jenkins’ method, but utilising PM7 to calculate the ∆H_f(g) values for the ions, is 112.6 kJ mol⁻¹, showing that the earlier result was in error due to computation of the ∆H_f(g) terms for the molecular ions.

For the 11 energetic co-crystals investigated, ∆H_L was calculated using two different approaches, the quicker PIXEL method, and the more computationally demanding DFT-D method (see Figure 3C and Supplementary Table 6), as no experimental data was available to benchmark the predictions against. From these results, it was readily apparent that the two methods provide comparable results (R² = 0.997, m = 0.964), meaning that both options are viable for co-crystals as part of a computational screening program. PIXEL does present some limitations, however, as it is not applicable to structures where the number of molecules in the crystallographic asymmetric unit exceeds 2; for these larger crystal lattices DFT-D is at present the only realistic solution.

Next, we offer an interesting comparison between the ∆H_L of the co-crystals (from Eq. 9) alongside the ∆H_L values for their respective co-formers (from Eqs. 5, 10). This is shown in Figure 4 (and Supplementary Tables 1, 6), from which it is apparent that the sum of the latter gives a reasonable approximation of the former. This is in agreement with Vener et al., who studied a number of co-crystals using DFT-D and hybrid-DFT functionals (Vener et al., 2014). A similar observation was also made by Day et al., who showed using DFT-D simulations that the total energy of 350 organic co-crystals was greater than a linear weighed sum of their single-component counterparts by the order of just ca. 8 kJ mol⁻¹ (per molecule). This applied to around 95% of their data set, showing that co-crystal formation is overwhelmingly driven by thermodynamics (Taylor and Day, 2018). Estimating ∆H_L for a co-crystal by this quick approximation may help guide the choice of co-formers for EM research, in order to maximize the amount of stored chemical energy. It also would be particularly useful to estimate ∆H_L for co-crystals such as CL-20/HMX, which has a 2:1 ratio of co-formers (Bolton et al., 2012) and thus falls outside the scope of PIXEL calculations. Accordingly, we estimate ∆H_L of CL-20/HMX to be ca. −430 kJ mol⁻¹ [based on ∆H_L (HMX) + 2 × ∆H_L (CL-20), see Supplementary Table 1]. Using DFT-D, the corresponding ∆H_L prediction is -460 kJ mol⁻¹. This results in a predicted ∆H_f(s) of ca. 570–600 kJ mol⁻¹. We note that Zhang et al. (2019) estimated a considerably higher value for ∆H_f(s) of 861.9 ± 18.6 kJ mol⁻¹; this was obtained indirectly using a calculated heat of reaction for the formation of solid CL-20/HMX from solid ε-CL-20 and β-HMX, that in turn used a thermochemical cycle based on measured enthalpies of dissolution of all species in acetonitrile, and literature ∆H_f(s) values for ε-CL-20 and β-HMX.

FIGURE 4

FIGURE 4. Comparison of calculated lattice enthalpies of co-crystals 68–78 by (A) PIXEL and (B) DFT-D, alongside constituent co-formers.

While Day et al. makes a strong case for co-crystallisation being overwhelmingly thermodynamically driven (Taylor and Day, 2018), a recent report by Perlovich suggests that the formation of around 30% of co-crystals are entropically driven (Perlovich, 2020). This work was based on constructing a dataset of 1947 co-crystals for which experimental sublimation energies were available, from which an algorithm to calculate the Gibbs sublimation energy was derived. Our test set of eleven EM co-crystals is too small to add any significant weight to this argument. While in general the predictions show that the lattice energy of the co-crystal does indeed exceed that of the sum of the co-formers (by up to ca. 70 kJ mol⁻¹, a similar total energy range reported by Day et al.), for some structures (68, 74, 76 and 78 by PIXEL, 68, 72 and 76 by DFT-D) the relationship does not hold, although we note that the energy shortfalls are typically very small. Improving the computational method to improve the accuracy of the comparison (i.e., zero point energy and entropy corrections) (Nyman and Day, 2015) will be considered more broadly in follow up studies. Such corrections could be obtained from calculating full phonon spectra, but these are time consuming to perform and therefore go against the computational screening philosophy which pervades this paper.

In summary, of the three methods tested here for calculating ∆H_f(g), PM7 performs favourably compared with the isodesmic reactions and atom equivalence methods, while offering application to the widest range of molecules. Its strength has also been demonstrated when calculating ∆H_f(s) for single component solids and salts, by the ESP and Jenkins methods, respectively, with predicted values showing excellent correlation with experimental ∆H_f(s) values. Finally, two different methods, PIXEL and DFT-D, were compared for calculating ∆H_L terms for co-crystals and found to give comparable results. Comparison of the calculated ∆H_L terms show that, to a first approximation, the lattice enthalpy of the co-crystal is the sum of the lattice enthalpies of the constituent co-formers. While this does not inform on the likelihood of success for the formation of new co-crystals, it does offer significant new insight in directing co-crystallisation studies to create new EMs with desired ∆H_f(s) values.

Local Force-Constant Calculations

Local vibrational-mode force constants were calculated for all CHNO-containing molecules in the test set, with a further ten EMs added to provide a wider and more comprehensive coverage of bond length values (31 CHNO molecules in total, Supplementary Table 7, including the EMs CL-20, RDX, HMX, HNB, NTO, TATB, FOX-7, PETN and nitroglycerin). From this data, a relationship between the bond lengths and force constants (bond strengths) of seven different covalent bonding environments can be drawn (see Figure 5).

FIGURE 5

FIGURE 5. Correlation between bond lengths and local vibrational mode force constants for seven covalent bonding environments.

The data shown here expands upon the relationship that has been established by Kraka et al. for C−C bonds occurring in both gas and solid state geometries.(Kraka et al., 2020). It also mirrors the trends shown by Byler et al. (1987) and Ladd et al. (1964) for C−N bonds and C−O bonds, respectively. The studies by Byler and Ladd analyzed mainly linear molecules in an effort to reduce the effect of coupling of normal vibrational modes. This limitation is elegantly side-stepped by local vibrational mode analysis through the recasting of the normal modes of vibration into the local modes through mass decoupling (Groom et al., 2016; Hehre et al., 1970) Byler et al. reported bond lengths that ranged from 1.122 Å (k = 20.17 mdyn Å⁻¹, H₃CC−NBCl₃) to 1.555 Å (k = 3.51 mdyn Å⁻¹, F₃C−NO). Their data sit comfortably on our C−N curve, highlighting the power of the local-mode analysis route for obtaining force constants for more complex molecules. Ladd et al. measured the vibrational frequency of the C−O bond stretch for different excited states of carbon monoxide to obtain the relationship between force constant and bond length over the range 1.088–1.396 Å. Their data agree quantitatively with the C–O bond curve shown in Figure 5 at long bond lengths, but consistently underestimates the force constants at shorter (<1.2 Å) distances, which suggests that measuring the force constants of only excited states of CO has skewed the relationship between bond length and force constant. McKean made extensive studies of the vibrational frequencies of isolated C–H bond stretches, a property which is directly comparable to the local mode of vibration (McKean, 1978; Konkoli et al., 1998; Konkoli and Cremer, 1998a; Larsson and Cremer, 1999), and therefore the force constants, and presented a relationship between C–H bond length and experimentally determined stretching frequencies which mirrors the trend shown here. Cremer et al. (2000) also calculated C–H force constants of adiabatic internal modes (an earlier name for local vibrational modes), their shortest C–H bond being 1.086 Å with a force constant of 5.58 mdyn Å⁻¹. This fits on our line presented in Figure 5, which is now further extended to 1.066 Å (and 6.46 mdyn Å⁻¹).

A strong correlation of force constant and bond length is observed for all bond types, and decay functions were used to fit trend lines (Supplementary Table 8) that returned R² values of ca. 0.99 for all bond types, with the exception of C−N (R² = 0.98) and N−H (R² = 0.92), although the latter corresponds to only four data points and is likely under-represented. Four bond types, C−C, C−N, C−O and N−N, encompass single to triple bond behaviour, and N−O bond types include single and double bonds, as evident from the clustering of data points. The wide range of bond lengths studied give rise to a corresponding wide range of calculated local-force constant values, which for the most part follow the sequence C−C > C−N ≈ C−O > N−N ≈ N−O, although the longer distances associated with single C−N bonds render these interactions on a par with single N−N and N−O bonds. The weakness of the long C–N bond fits with expectation: the rupture of R–NO₂ bonds has been shown to be a critical step in the decomposition of energetic materials (Sharma and Owens, 1979; Sharma et al., 1982; Owens, 1996). The data shown in Figure 5 are testament to the great potential offered by the local-mode analysis route: it is a quick and straightforward method to compare and contrast the bond strengths of all bonding interactions within a molecule, from which the weakest bond can be unambiguously identified. In essence, the curves shown in Figure 5 allow a direct mapping between bond length and bond strength to be obtained. For the case of CL-20, a cage nitramine structure comprising 5, 6 and 7 membered C−N rings (compound 21 in Figure 1), it has been suggested that the C−N bonds forming the cage can also act as the “trigger linkage,” (Sun et al., 2021) in addition to the generally weaker N−NO₂ bonds. This is supported here, with the C−N bonds in the more strained 5-membered rings having force constants similar to the stronger N−NO₂ bonds (3.337 vs. 3.151 mdyn Å⁻¹, respectively). The less strained 6-membered ring contains stronger C−N bonds, with force constants ranging from 4.088–4.486 mdyn Å⁻¹, a considerable increase in strength from the C−N bonds present in the 5-membered rings. The relationship between bond length and force constant presented here therefore has the potential to be applied to molecular design, as the weakest bond in a molecule can be tuned by its surrounding molecular environment.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

IC performed the calculations. IC and CM wrote the manuscript. All authors analysed the data. AM, CP, and CM directed the research.

Funding

This material is based upon work supported by the Air Force Office of Scientific Research under award number FA8655-20-1-7000. We are grateful for computational support from the United Kingdom Materials and Molecular Modelling Hub, which is partially funded by EPSRC (EP/P020194 and EP/T022213), for which access was obtained via the UKCP consortium and funded by EPSRC grant ref EP/P022561/1.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling Editor declared a past co-authorship with the authors (AM, CP, CM)

Acknowledgments

We gratefully acknowledge Prof. E. Kraka (Department of Chemistry, Southern Methodist University, Texas) for early access to their local mode analysis software, LMODEA. IC thanks the University of Edinburgh, School of Chemistry for the award of a PhD studentship. We thank the EaSTCHEM Research Computing Facility and the Edinburgh Computing Data Facility, for access to software and hardware resources, respectively.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem.2021.726357/full#supplementary-material

References

Axthammer, Q. J., Krumm, B., and Klapötke, T. M. (2015). Pentaerythritol-based Energetic Materials Related to PETN. Eur. J. Org. Chem. 2015, 723–729. doi:10.1002/ejoc.201403265

CrossRef Full Text | Google Scholar

Benson, S. W., Cruickshank, F. R., Golden, D. M., Haugen, G. R., O'Neal, H. E., Rodgers, A. S., et al. (1969). Additivity Rules for the Estimation of Thermochemical Properties. Chem. Rev. 69, 279–324. doi:10.1021/cr60259a002

CrossRef Full Text | Google Scholar

Bolton, O., Simke, L. R., Pagoria, P. F., and Matzger, A. J. (2012). High Power Explosive with Good Sensitivity: A 2:1 Cocrystal of CL-20:HMX. Cryst. Growth Des. 12, 4311–4314. doi:10.1021/cg3010882

CrossRef Full Text | Google Scholar

Bozkuş, S. I., Hope, K. S., Yüksel, B., Atҫeken, N., Nazır, H., Atakol, O., et al. (2019). Characterization and Properties of a Novel Energetic Co-crystal Formed between 2,4,6-Trinitrophenol and 9-Bromoanthracene. J. Mol. Struct. 1192, 145–153. doi:10.1016/j.molstruc.2019.04.109

CrossRef Full Text | Google Scholar

Brandenburg, J. G., Maas, T., and Grimme, S. (2015). Benchmarking DFT and Semiempirical Methods on Structures and Lattice Energies for Ten Ice Polymorphs. J. Chem. Phys. 142, 124104–124111. doi:10.1063/1.4916070

CrossRef Full Text | Google Scholar

Byler, D. M., Susi, H., and Damert, W. C. (1987). Relation between Force Constant and Bond Length for Carbon-Nitrogen Bonds. Spectrochimica Acta A: Mol. Spectrosc. 43, 861–863. doi:10.1016/0584-8539(87)80231-3

CrossRef Full Text | Google Scholar

Byrd, E. F. C., and Rice, B. M. (2006). Improved Prediction of Heats of Formation of Energetic Materials Using Quantum Mechanical Calculations. J. Phys. Chem. A. 110, 1005–1013. doi:10.1021/jp0536192

CrossRef Full Text | Google Scholar

Chan, B., Collins, E., and Raghavachari, K. (2020). Applications of Isodesmic‐type Reactions for Computational Thermochemistry. Wires Comput. Mol. Sci. 11, 1–18. doi:10.1002/wcms.1501

CrossRef Full Text | Google Scholar

Clark, S. J., Segall, M. D., Pickard, C. J., Hasnip, P. J., Probert, M. I. J., Refson, K., et al. (2005). First Principles Methods Using CASTEP. Z. fuer Krist 220, 567–570. doi:10.1524/zkri.220.5.567.65075

CrossRef Full Text | Google Scholar

Cramer, C. (2004). Essentials of Computational Chemistry Theories and Models. 2nd ed. Chichester: Wiley.

Cremer, D., Wu, A., Larsson, A., and Kraka, E. (2000). Some Thoughts about Bond Energies, Bond Lengths, and Force Constants. J. Mol. Model. 6, 396–412. doi:10.1007/PL00010739

CrossRef Full Text | Google Scholar

Cumming, A. S. (2017). Environmental Aspects of Energetic Materials Use and Disposal. Chem. Rocket Propulsion, 727–741. doi:10.1007/978-3-319-27748-6

CrossRef Full Text | Google Scholar

Curtiss, L. A., Raghavachari, K., Redfern, P. C., and Pople, J. A. (1997). Assessment of Gaussian-2 and Density Functional Theories for the Computation of Enthalpies of Formation. J. Chem. Phys. 106, 1063–1079. doi:10.1063/1.473182

CrossRef Full Text | Google Scholar

Elioff, M. S., Hoy, J., and Bumpus, J. A. (2016). Calculating Heat of Formation Values of Energetic Compounds: A Comparative Study. Adv. Phys. Chem., 2016, 1. doi:10.1155/2016/5082084

CrossRef Full Text | Google Scholar

Fan, J.-Y., Zheng, Z.-Y., Su, Y., and Zhao, J.-J. (2017). Assessment of Dispersion Correction Methods within Density Functional Theory for Energetic Materials. Mol. Simulation 43, 568–574. doi:10.1080/08927022.2017.1293258

CrossRef Full Text | Google Scholar

Feng, S., and Li, T. (2006). Predicting Lattice Energy of Organic Crystals by Density Functional Theory with Empirically Corrected Dispersion Energy. J. Chem. Theor. Comput. 2, 149–156. doi:10.1021/ct050189a

PubMed Abstract | CrossRef Full Text | Google Scholar

Fomin, V. N., Gogol, D. B., Rozhkovoy, I. E., and Ponomarev, D. L. (2017). Quantum Chemical and Thermodynamic Calculations of Fulvic and Humic Copper Complexes in Reactions of Malachite and Azurite Formation. Appl. Geochem. 79, 9–16. doi:10.1016/j.apgeochem.2017.02.002

CrossRef Full Text | Google Scholar

Fried, L., and Souers, P. (1994). CHEETAH: A Next Generation Thermochemical Code. doi:10.2172/95184

CrossRef Full Text | Google Scholar

Frisch, M. J., Trucks, G. W., Schlegel, H. B., Scuseria, G. E., Robb, M. A., Cheeseman, J. R., et al. (2016). Gaussian 16, Revision A.03.

Google Scholar

Gao, H., Ye, C., Piekarski, C. M., and Shreeve, J. n. M. (2007). Computational Characterization of Energetic Salts. J. Phys. Chem. C 111, 10718–10731. doi:10.1021/jp070702b

CrossRef Full Text | Google Scholar

Gavezzotti, A. (2005). Calculation of Lattice Energies of Organic Crystals: The PIXEL Integration Method in Comparison with More Traditional Methods. Z. Krist 220, 499–510. doi:10.1524/zkri.220.5.499.65063

CrossRef Full Text | Google Scholar

Gavezzotti, A. (2011). Efficient Computer Modeling of Organic Materials. The Atom-Atom, Coulomb-London-Pauli (AA-CLP) Model for Intermolecular Electrostatic-Polarization, Dispersion and Repulsion Energies. New J. Chem. 35, 1360–1368. doi:10.1039/c0nj00982b

CrossRef Full Text | Google Scholar

Groom, C. R., Bruno, I. J., Lightfoot, M. P., and Ward, S. C. (2016). The Cambridge Structural Database. Acta Crystallogr. Sect B 72, 171–179. doi:10.1107/S2052520616003954

CrossRef Full Text | Google Scholar

Hehre, W. J., Ditchfield, R., Radom, L., and Pople, J. A. (1970). Molecular Orbital Theory of the Electronic Structure of Organic Compounds. V. Molecular Theory of Bond Separation. J. Am. Chem. Soc. 92, 4796–4801. doi:10.1021/ja00719a006

CrossRef Full Text | Google Scholar

Jenkins, H. D. B., Roobottom, H. K., Passmore, J., and Glasser, L. (1999). Relationships Among Ionic Lattice Energies, Molecular (Formula Unit) Volumes, and Thermochemical Radii. Inorg. Chem. 38, 3609–3620. doi:10.1021/ic9812961

PubMed Abstract | CrossRef Full Text | Google Scholar

Jenkins, H. D. B., Tudela, D., and Glasser, L. (2002). Lattice Potential Energy Estimation for Complex Ionic Salts from Density Measurements. Inorg. Chem. 41, 2364–2367. doi:10.1021/ic011216k

PubMed Abstract | CrossRef Full Text | Google Scholar

Kennedy, S. R., and Pulham, C. R. (2018). “Chapter 6. Co-crystallization of Energetic Materials,” in In Co-crystals: Preparation, Characterization and Applications Edited. Editors B. Aakeroy, and A. Sinha (London, United Kingdom: The Royal Society of Chemistry), 231–266. doi:10.1039/9781788012874

CrossRef Full Text | Google Scholar

Konkoli, Z., and Cremer, D. (1998a). A New Way of Analyzing Vibrational Spectra. III. Characterization of normal Vibrational Modes in Terms of Internal Vibrational Modes. Int. J. Quant. Chem. 67, 29–40. doi:10.1002/(sici)1097-461x(1998)67:1<29::aid-qua3>3.0.co;2-0

CrossRef Full Text | Google Scholar

Konkoli, Z., and Cremer, D. (1998b). A New Way of Analyzing Vibrational Spectra. I. Derivation of Adiabatic Internal Modes. Int. J. Quant. Chem. 67, 1–9. doi:10.1002/(sici)1097-461x(1998)67:1<1::aid-qua1>3.0.co;2-z

CrossRef Full Text | Google Scholar

Konkoli, Z., Larsson, J. A., and Cremer, D. (1998). A New Way of Analyzing Vibrational Spectra. IV. Application and Testing of Adiabatic Modes within the Concept of the Characterization of normal Modes. Int. J. Quant. Chem. 67, 41–55. doi:10.1002/(sici)1097-461x(1998)67:1<41::aid-qua4>3.0.co;2-z

CrossRef Full Text | Google Scholar

Kraka, E., and Cremer, D. (2019). Dieter Cremer's Contribution to the Field of Theoretical Chemistry. Int. J. Quan. Chem. 119, e25849–28. doi:10.1002/qua.25849

CrossRef Full Text | Google Scholar

Kraka, E., Zou, W., and Tao, Y. (2020). Decoding Chemical Information from Vibrational Spectroscopy Data: Local Vibrational Mode Theory. Wires Comput. Mol. Sci. 10, 1–34. doi:10.1002/wcms.1480

CrossRef Full Text | Google Scholar

Ladd, J. A., Orville-Thomas, W. J., and Cox, B. C. (1964). Molecular Parameters and Bond Structure-III. Carbon-Oxygen Bonds. Spectrochimica Acta 20, 1771–1780. doi:10.1016/0371-1951(64)80180-6

CrossRef Full Text | Google Scholar

Larsson, J. A., and Cremer, D. (1999). Theoretical Verification and Extension of the McKean Relationship between Bond Lengths and Stretching Frequencies. J. Mol. Struct. 485-486, 385–407. doi:10.1016/S0022-2860(99)00093-9

CrossRef Full Text | Google Scholar

Lu, T., and Chen, F. (2012). Multiwfn: A Multifunctional Wavefunction Analyzer. J. Comput. Chem. 33, 580–592. doi:10.1002/jcc.22885

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, Q., Jiang, T., Chi, Y., Chen, Y., Wang, J., Huang, J., et al. (2017). A Novel Multi-Nitrogen 2,4,6,8,10,12-Hexanitrohexaazaisowurtzitane-Based Energetic Co-crystal with 1-Methyl-3,4,5-Trinitropyrazole as a Donor: Experimental and Theoretical Investigations of Intermolecular Interactions. New J. Chem. 41, 4165–4172. doi:10.1039/c6nj03976f

CrossRef Full Text | Google Scholar

Macrae, C. F., Bruno, I. J., Chisholm, J. A., Edgington, P. R., Mccabe, P., Pidcock, E., et al. (2008). Mercury CSD 2.0- New Features for the Visualization and Investigation of crystal Structures. J. Appl. Cryst. 41, 466–470. doi:10.1107/S0021889807067908

CrossRef Full Text | Google Scholar

Maloney, A. G. P., Wood, P. A., and Parsons, S. (2015). Intermolecular Interaction Energies in Transition Metal Coordination Compounds. CrystEngComm 17, 9300–9310. doi:10.1039/c5ce01522g

CrossRef Full Text | Google Scholar

McKean, D. C. (1978). Individual CH bond strengths in simple organic compounds: Effects of conformation and substitution. Chem Soc Rev. 7, 399–422. doi:10.1039/CS9780700399

CrossRef Full Text | Google Scholar

Michalchuk, A. A. L., Fincham, P. T., Portius, P., Pulham, C. R., and Morrison, C. A. (2018a). A Pathway to the Athermal Impact Initiation of Energetic Azides. J. Phys. Chem. C 122, 19395–19408. doi:10.1021/acs.jpcc.8b05285

CrossRef Full Text | Google Scholar

Michalchuk, A. A. L., Rudić, S., Pulham, C. R., and Morrison, C. A. (2018b). Vibrationally Induced Metallisation of the Energetic Azide α-NaN3. Phys. Chem. Chem. Phys. 20, 29061–29069. doi:10.1039/c8cp06161k

PubMed Abstract | CrossRef Full Text | Google Scholar

Michalchuk, A. A. L., Hemingway, J., and Morrison, C. A. (2021). Predicting the Impact Sensitivities of Energetic Materials through Zone-center Phonon Up-Pumping. J. Chem. Phys. 154, 064105–064111. doi:10.1063/5.0036927

CrossRef Full Text | Google Scholar

Michalchuk, A. A. L., Trestman, M., Rudić, S., Portius, P., Fincham, P. T., Pulham, C. R., et al. (2019). Predicting the Reactivity of Energetic Materials: an Ab Initio Multi-Phonon Approach. J. Mater. Chem. A. 7, 19539–19553. doi:10.1039/c9ta06209b

CrossRef Full Text | Google Scholar

Morrison, C. A., and Siddick, M. M. (2003). Determining the Strengths of Hydrogen Bonds in Solid-State Ammonia and Urea: Insight from Periodic DFT Calculations. Chem. Eur. J. 9, 628–634. doi:10.1002/chem.200390067

PubMed Abstract | CrossRef Full Text | Google Scholar

Nyman, J., and Day, G. M. (2015). Static and Lattice Vibrational Energy Differences between Polymorphs. CrystEngComm 17, 5154–5165. doi:10.1039/c5ce00045a

CrossRef Full Text | Google Scholar

Owens, F. J. (1996). Calculation of Energy Barriers for Bond Rupture in Some Energetic Molecules. J. Mol. Struct. THEOCHEM 370, 11–16. doi:10.1016/S0166-1280(96)04673-8

CrossRef Full Text | Google Scholar

Perdew, J. P., Burke, K., and Ernzerhof, M. (1996). Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 77, 3865–3868. doi:10.1103/PhysRevLett.77.3865

PubMed Abstract | CrossRef Full Text | Google Scholar

Perlovich, G. L. (2020). Formation Thermodynamics of Two-Component Molecular Crystals: Polymorphism, Stoichiometry, and Impact of Enantiomers. Cryst. Growth Des. 20, 5526–5537. doi:10.1021/acs.cgd.0c00695

CrossRef Full Text | Google Scholar

Piercey, D. G., Chavez, D. E., Scott, B. L., Imler, G. H., and Parrish, D. A. (2016). An Energetic Triazolo-1,2,4-Triazine and its N-Oxide. Angew. Chem. Int. Ed. 55, 15315–15318. doi:10.1002/anie.201608723

CrossRef Full Text | Google Scholar

Politzer, P., Lane, P., and Concha, M. C. (2003). Computational Approaches to Heats of Formation. Elsevier B.V 12, 247–277. doi:10.1016/s1380-7323(03)80011-0

CrossRef Full Text | Google Scholar

Politzer, P., Murray, J. S., Edward Grice, M., Desalvo, M., and Miller, E. (1997). Calculation of Heats of Sublimation and Solid Phase Heats of Formation. Mol. Phys. 91, 923–928. doi:10.1080/002689797171030

CrossRef Full Text | Google Scholar

Politzer, P., and Murray, J. S. (1998). Relationships between Lattice Energies and Surface Electrostatic Potentials and Areas of Anions. J. Phys. Chem. A. 102, 1018–1020. doi:10.1021/jp972885f

CrossRef Full Text | Google Scholar

Politzer, P., and Murray, J. S. (2002). The Fundamental Nature and Role of the Electrostatic Potential in Atoms and Molecules. Theor. Chem. Acc. Theor. Comput. Model. (Theoretica Chim. Acta) 108, 134–142. doi:10.1007/s00214-002-0363-9

CrossRef Full Text | Google Scholar

Politzer, P., Seminario, J. M., and Concha, M. C. (1998). Statitical Analysis of the Molecular Surfae Electrostatic Potential: an Approach to Describing Noncovalent Interactions in Condensed Phases. J. Mol. Struct. 427, 123–129. doi:10.1016/S0166-1280(97)00162-0

CrossRef Full Text | Google Scholar

Ponomarev, D. A., and Takhistov, V. V. (1997). What Are Isodesmic Reactions? J. Chem. Educ. 74, 201–203. doi:10.1021/ed074p201

CrossRef Full Text | Google Scholar

Reeves, M. G., Wood, P. A., and Parsons, S. (2020). MrPIXEL: Automated Execution of Pixel Calculations via theMercuryinterface. J. Appl. Cryst. 53, 1154–1162. doi:10.1107/S1600576720008444

PubMed Abstract | CrossRef Full Text | Google Scholar

Sabatini, J., and Oyler, K. (2015). Recent Advances in the Synthesis of High Explosive Materials. Crystals 6, 5–22. doi:10.3390/cryst6010005

CrossRef Full Text | Google Scholar

Sharma, J., Garrett, W. L., Owens, F. J., and Vogel, V. L. (1982). X-ray Photoelectron Study of Electronic Structure, Ultraviolet, and Isothermal Decomposition of 1,3,5-Triamino-2,4,6-Trinitrobenzene. J. Phys. Chem. 86, 1657–1661. doi:10.1021/j100206a034

CrossRef Full Text | Google Scholar

Sharma, J., and Owens, F. J. (1979). XPS Study of UV and Shock Decomposed Triamino-Trinitrobenzene. Chem. Phys. Lett. 61, 280–282. doi:10.1016/0009-2614(79)80644-2

CrossRef Full Text | Google Scholar

Sinditskii, V. P., Filatov, S. A., Kolesov, V. I., Kapranov, K. O., Asachenko, A. F., Nechaev, M. S., et al. (2015). Combustion Behavior and Physico-Chemical Properties of Dihydroxylammonium 5,5′-Bistetrazole-1,1′-Diolate (TKX-50). Thermochim. Acta 614, 85–92. doi:10.1016/j.tca.2015.06.019

CrossRef Full Text | Google Scholar

Singh, H. J., Gupta, S., and Sengupta, S. K. (2014a). Computational Studies on Nitramino Derivatives of 1-Amino-1,2-Azaboriridine as High Energetic Material. RSC Adv. 4, 40534–40541. doi:10.1039/c4ra06285j

CrossRef Full Text | Google Scholar

Singh, H. J., Upadhyay, M. K., and Sengupta, S. K. (2014b). Theoretical Studies on Benzo[1,2,4]triazine-Based High-Energy Materials. J. Mol. Model. 20, 1–10. doi:10.1007/s00894-014-2205-9

CrossRef Full Text | Google Scholar

Stewart, J. J. P. (2013). Optimization of Parameters for Semiempirical Methods VI: More Modifications to the NDDO Approximations and Re-optimization of Parameters. J. Mol. Model. 19, 1–32. doi:10.1007/s00894-012-1667-x

CrossRef Full Text | Google Scholar

Suceska, M. (2004). Calculation of Detonation Parameters by EXPLO5 Computer Program. Mater. Sci. Forum 465–466, 325–330. doi:10.4028/www.scientific.net/msf.465-466.325

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, K.-b., Zhang, S.-h., Ren, F.-d., Hao, Y.-P., and Ba, S.-h. (2021). Theoretical Prediction of the Trigger Linkage, Cage Strain, and Explosive Sensitivity of CL-20 in the External Electric fields. J. Mol. Model. 27, 1–20. doi:10.1007/s00894-020-04634-8

CrossRef Full Text | Google Scholar

Taylor, C. R., and Day, G. M. (2018). Evaluating the Energetic Driving Force for Cocrystal Formation. Cryst. Growth Des. 18, 892–904. doi:10.1021/acs.cgd.7b01375

PubMed Abstract | CrossRef Full Text | Google Scholar

Vatani, A., Mehrpooya, M., and Gharagheizi, F. (2007). Prediction of Standard Enthalpy of Formation by a QSPR Model. Int. J. Mol. Sci. 8, 407–432. doi:10.3390/i8050407

CrossRef Full Text | Google Scholar

Vener, M. V., Levina, E. O., Koloskov, O. A., Rykounov, A. A., Voronin, A. P., and Tsirelson, V. G. (2014). Evaluation of the Lattice Energy of the Two-Component Molecular Crystals Using Solid-State Density Functional Theory. Cryst. Growth Des. 14, 4997–5003. doi:10.1021/cg5005243

CrossRef Full Text | Google Scholar

Wan, Z., Wang, Q. D., and Liang, J. (2020). Accurate Prediction of Standard Enthalpy of Formation Based on Semiempirical Quantum Chemistry Methods with Artificial Neural Network and Molecular Descriptors. Int. J. Quan. Chem. 121, 1–16. doi:10.1002/qua.26441

CrossRef Full Text | Google Scholar

Zhang, S., Zhang, J., Kou, K., Jia, Q., Xu, Y., Liu, N., et al. (2019). Standard Enthalpy of Formation, Thermal Behavior, and Specific Heat Capacity of 2HNIW·HMX Co-crystals. J. Chem. Eng. Data 64, 42–50. doi:10.1021/acs.jced.8b00454

CrossRef Full Text | Google Scholar

Zhang, Z.-B., Li, T., Yin, L., Yin, X., and Zhang, J.-G. (2016). A Novel Insensitive Cocrystal Explosive BTO/ATZ: Preparation and Performance. RSC Adv. 6, 76075–76083. doi:10.1039/c6ra14510h

CrossRef Full Text | Google Scholar

Zou, W., Kalescky, R., Kraka, E., and Cremer, D. (2012). Relating normal Vibrational Modes to Local Vibrational Modes with the Help of an Adiabatic Connection Scheme. J. Chem. Phys. 137, 084114. doi:10.1063/1.4747339

CrossRef Full Text | Google Scholar