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. Δ H R = Δ E 0 + Δ ZPE + Δ H T + Δ 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. Δ H R = ∑ Δ H f ( g ) p r o d u c t s   + ∑ Δ H f ( g ) r e a c t a n t s   (2)

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). Δ H f ( g ) = E − ∑ n j e j   (3)

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.

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. Δ H f ( s )   = Δ H f ( g )   − Δ H s u b (4) Δ H s u b = a ( S A ) 2 + b σ t o t 2 ν   + c (5) Where a ,   b and c are semi-empirically deduced fitting parameters proposed by Byrd et al. (Byrd and Rice, 2006) S A is the surface area of the 0.001 electron.bohr⁻³ isosurface of the electrostatic potential of the molecule, σ t o t 2   is the measure of variablity of electronic potential on the surface, and ν 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).

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. Δ H f ( s )   =   Δ H f ( g ) ( c a t i o n ) +   Δ H f ( g ) ( a n i o n ) −   Δ H L (6)

Here, ∆H_L is expressed by Eq. 7, as proposed by Jenkins et al. (2002). Δ H L = U p o t + [ p ( n m / 2 − 2 ) + q ( n x   / 2 − 2 ) ] R T (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 p o t denotes the lattice potential energy and in turn is defined by Eq. 8. U p o t = α ( V m ) 1 3 + β (8) Where V m   denotes the molecular volume (V_cell/Z), in nm³, and the remaining coefficients α and   β are fitting terms provided by Jenkins et al. and which vary depending on the charge ratio of the salt.

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. Δ H L = E c e l l Z − E c o − f o r m e r 1 − E c o − f o m e r 2 (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. Δ H L =   − Δ H s u b − 2 R T (10)

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.

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.

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.

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.

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.

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.

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

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.