Calculating Heat of Formation Values of Energetic Compounds : A Comparative Study

Heat of formation is one of several important parameters used to assess the performance of energetic compounds. We evaluated the ability of six different methods to accurately calculate gas-phase heat of formation (ΔfH o 298,g) values for a test set of 45 nitrogencontaining energetic compounds. Density functional theory coupled with the use of isodesmic or other balanced equations yielded calculated results in which 82% (37 of 45) of the ΔfH o 298,g values were within ±2.0 kcal/mol of the most recently recommended experimental/reference values available. This was compared to a procedure using density functional theory (DFT) coupled with an atom and group contribution method in which 51% (23 of 45) of the ΔfH o 298,g values were within ±2.0 kcal/mol of these values.The T1 procedure andBenson’s group additivitymethod yielded results inwhich 51% (23 of 45) and 64% (23 of 36) of theΔfH o 298,g values, respectively, were within ±2.0 kcal/mol of these values. We also compared two relatively new semiempirical approaches (PM7 and RM1) with regard to their ability to accurately calculate ΔfH o 298,g. Although semiempirical methods continue to improve, they were found to be less accurate than the other approaches for the test set used in this investigation.


Introduction
There is substantial interest in the discovery and development of new energetic compounds which include high explosives and propellants and there is special interest in the discovery and development of high energy density materials (HEDMs). Collectively, compounds that have densities greater than or equal to 1.9 g/cm 3 , detonation velocities greater than or equal to 9.0 km/s, and detonation pressures greater than or equal to 40.0 GPa are known as HEDMs [1]. The usefulness of HEDMs in a variety of processes has been understood for at least a century [2]. Two important considerations with regard to HEDMs are that they can be dangerous and costly to synthesize. Moreover, increasing environmental concerns call for more effective ways of predicting performance of HEDMs [3]. Manufacturers must be able to determine the amount of energy that can be stored in a molecule while maintaining an acceptable level of stability [4].
The condensed-phase heat of formation (Δ 298,s or Δ 298,l ) is one of several important parameters used to assess the performance of energetic materials [5]. In practice, most theoretical approaches calculate gas-phase heat of formation (Δ 298,g ) values. Solid-or liquid-phase values are then calculated by subtracting heat of sublimation (Δ 298,sub ) or heat of vaporization (Δ 298,vap ) values, respectively, from the gas-phase value.
Our interest concerns the use of relatively rapid methods for the calculation/prediction of accurate gas-phase heat of formation values. During the past several years two relatively rapid computational procedures have been described for use in determining Δ 298,g values [6][7][8]. One procedure, developed by Byrd and Rice [6,7], uses density functional theory (DFT) coupled with an atom and group contribution method to determine Δ 298,g for energetic compounds. Another procedure, developed by Ohlinger et al. [8], was designed to be applicable to a broader range of compounds. The latter relies on a novel multilevel computational approach known as T1 that approaches the accuracy of the G3MP2 method [9] while simultaneously reducing the computation time by 2-3 orders of magnitude [8]. In the investigations reported here, 2 Advances in Physical Chemistry we have compared the effectiveness of these two procedures with regard to their ability to accurately calculate Δ 298,g values of nitrogen-containing energetic compounds. Operating under the assumption that newer procedures are not necessarily more reliable than more established methods, we also assessed and compared the ability of two other wellestablished procedures used to calculate Δ 298,g . Specifically, we assessed and compared results using the group additivity method developed by Benson and his colleagues [10,11] as well as density functional theory coupled with the use of isodesmic, isogyric, or other balanced equations [12]. We also used and compared two relatively new semiempirical models (PM7 and RM1) [13,14] to calculate Δ 298,g .

Methods
The test set used in this investigation consisted of the gasphase heat of formation (Δ 298,g ) values for forty-five nitrogen-containing energetic compounds studied by Byrd and Rice [6,7]. Δ 298,g values were calculated using the T1 procedure which is based on the G3MP2 [9] multilevel ab initio model chemistry as implemented in the Spartan'10 and 14 suite of programs. Semiempirical theory using the RM1 and PM7 model chemistries [13,14] were used to calculate Δ 298,g as implemented in Spartan 14 (Wavefunction Inc., Irvine, CA) and MOPAC2012 (http://openmopac.net/), respectively.
Gas-phase heat of formation values were also calculated using density functional theory as implemented in Spartan 10 and 14 coupled with the use of isodesmic, isogyric and other balanced equations. An informative account of the details concerning the calculation of Δ 298,g using this approach has been presented by Cramer [12]. Most of the equations used were isodesmic equations. For the studies reported herein total electronic energies ( Tot ), zero point energies (ZPE), and thermal correction ( ) values for compounds in all of the equations were calculated using the M06 2X functional [15,16] and the 6-31G * basis set. The 6-31G * basis set is a medium size basis set. It was selected for use so that excessive computation time could be avoided. In Spartan, it is necessary that the "compute I.R." box be selected in order to make these calculations. When this is done, the calculated 298 value required for use in isodesmic, isogyric, or other balanced equations is computed and appears in the Thermodynamics Window. This window is accessed by the following path: Display > Properties > Thermodynamics. Experimental values for Δ 298,g required for compounds used in these studies were acquired from the NIST WebBook [17], Pedley's Thermochemical Data of Organic Compounds, 2nd ed. [18], Cox and Pilcher's, Thermochemistry of Organic and Organometallic Compounds, [19] and publications by Dorofeeva and her colleagues [20][21][22][23][24].
Gas-phase heat of formation values were calculated for twenty-five compounds using the Group Additivity approach as implemented in the NIST WebBook (http://webbook.nist .gov/chemistry/) applet [17]. This applet uses Benson Group Additivity values [10,11]. Eight compounds in our test set contained an azido group. Unfortunately, the azido group value is not among those included in the database for the NIST WebBook applet. To address this problem we first calculated an enthalpy group value for an azido group attached to a carbon atom. The Δ 298,g of methyl azide has been determined experimentally and by high level computation. We used the value of 71.2 kcal/mol (±1.0 kcal/mol) recommended by Dorofeeva et al. [21]. Subtracting the methyl group value (−10.2 kcal/mol) from the heat of formation of methyl azide gives a value of 81.4 kcal/mol for the azido group. The group additivity approach was then used to calculate Δ 298,g values of azido-containing compounds as follows. First the NIST applet was used to calculate the Δ 298,g value for the amino analog of the azido compound. Then, the azido group value of 81.4 kcal/mol was used in the calculation instead of the amino group value. Using this simple modification of the group additivity approach Δ 298,g values were calculated for six of the eight compounds in our test set that contained an azido group. Calculation of Δ 298,g for the remaining compounds that contained an azido group is described below.
The ASTM CHETAH 8.0 software program [25] was used to calculate Δ 298,g values for a few other compounds in our test set whose Δ 298,g values could not be calculated using the NIST applet.
The group additivity difference method as described by Cohen and Benson [26] was used to calculate the Δ 298,g value for nitromethane. For nitromethane, it should be mentioned that experimental values of −17.8 and −19.3 kcal/mol have been reported and used extensively. Interestingly, there is considerable theoretical as well as experimental support for both values [20,[22][23][24][27][28][29]. In the present investigation, the group additivity value for nitromethane (−18.0 kcal/mol) was calculated by subtracting the group value (−6.6 kcal/mol) for a methylene (-CH 2 -) group bound to a nitrogen atom from the calculated Δ 298,g value of −24.6 kcal/mol for nitroethane. We note that this value is 1.5 kcal/mol greater than that previously reported by Bumpus and Willoughby [27]. It should also be noted that recent experimental results [24] support a slightly greater Δ 298,g value of −17.09 kcal/mol for nitromethane and it is this value that was ultimately included in the alternative experimental reference set.
The calculation of Δ 298,g for 1-azido-1,1-dinitroethane was accomplished using the group values for the methyl group and azido group coupled with the group value for a carbon attached to an azido group and two nitro groups. This later value was assumed to be equivalent to the group value (−9.9 kcal/mol) for a carbon attached to a hydrogen and two nitro groups.
For azidotrinitromethane, the group value for the trinitromethyl group (−1.45 kcal/mol) listed in CHETAH 8.0 was added to the azido group value of 81.4 kcal/mol giving a calculated Δ 298,g value of 79.95 kcal/mol (= 80.0 kcal/mol) for azidotrinitromethane.
All calculated Δ 298,g values were compared with experimental values and with values published by Byrd and Rice [6,7] who calculated Δ 298,g values using density Advances in Physical Chemistry 3 functional theory coupled with an atom and group equivalent method. Calculated Δ 298,g values were also compared to an alternative experimental reference set in which eighteen presumably more accurate Δ 298,g values recommended by Dorofeeva and colleagues [20][21][22][23][24] were substituted for corresponding values in the original set of experimental values. Also substituted in this alternative experimental reference set was the Δ 298,g value for 2,4,6-trinitrotoluene listed in Cox and Pilcher's "Thermochemistry of Organic and Organometallic Compounds" [19]. Finally, it should be noted that during the course of this investigation it appeared that the experimental Δ 298,g value for dinitromethylbenzene might be several kcal/mol too low. To address this possibility we used an approach similar to that used successfully by Dorofeeva et al. [21][22][23] to deal with such issues. Specifically we used density functional theory (EDF2 [30] functional and the 6-311++G (2df,2p) basis set) coupled with the use of ten isodesmic equations to calculate a Δ 298,g reference value of 14.2 ± 0.3 for dinitromethylbenzene. This value was also included in the alternative reference set.
It should be mentioned that many of the computed RM1 and PM7 Δ 298,g values for compounds appearing in Table 1 have been previously reported [13,14]. Similarly most of the T1 Δ 298,g values appear in the Spartan Spectra and Properties Database (SSPD) data base that accompanies Spartan 14. In general, our results are quite consistent with these previously reported values. In some instances, however, slightly lower Δ 298,g values were calculated indicating the existence of lower-energy equilibrium geometries than those used in previous investigations.

Results and Discussion
The ultimate goal of any computational strategy designed to estimate values for Δ 298,g is to reduce uncertainties to within ±1.0 or 2.0 kcal/mol (1 kcal/mol = 4.184 kJ/mol), which is taken to be experimental error. One reason for this goal is that small errors in Δ 298,g become magnified when propagated to obtain other values which depend on Δ 298,g . The standard Gibbs energy change of a reaction, for example, is a function of the reaction enthalpy change as follows: Reaction enthalpies can be calculated in accord with Hess's law by taking the difference of the standard heat of formation values of the products and reactants of the reaction. The equilibrium constant ( eq ) of a reaction and Δ rxn ∘ are related by the expression Since the equilibrium constant depends exponentially upon Δ rxn ∘ , a small error in Δ rxn can result in dramatically different values for eq when it is the calculated parameter of interest.
For energetic compounds, calculation of equilibrium constants using heat of formation values is typically not a concern since most reactions of interest are highly exergonic with equilibrium constants that indicate product formation is highly favored. It must be noted, however, that several parameters (e.g., heat of detonation, heat of explosion, power index, explosive velocity, and explosive pressure) used to characterize high explosives require the condensed-phase heat of formation value. This is somewhat problematic for theoretical compounds since calculations are performed on single molecules, which are, by definition, gas-phase molecules. To calculate solid-phase heat of formation values, heat of sublimation (or vaporization) values are calculated by one of several computational methods available and these values are then subtracted from the gas-phase values to calculate the condensed-phase values. This represents still another source of uncertainty in subsequent calculations. It is therefore important to be able to identify and use procedures that provide the most accurate heat of formation values possible.
Our initial objective was to compare the abilities of the T1 multilevel ab initio approach developed by Ohlinger et al. [8] with the atom and group equivalent DFT approach described by Byrd and Rice [6,7]. To make this comparison we selected the forty-five nitrogen-containing compounds studied by Byrd and Rice [6,7] which we designated as our test set and acquired T1 Δ 298,g values for these compounds. These results and the results of all comparisons made in this investigation are presented in Tables 1, 2, and 3 and in Figure 1.
When making comparisons such as those under consideration here, it is necessary to use the most accurate experimental Δ 298,g reference values available. This is problematic for no fewer than twenty of the experimental Δ 298,g values used in this investigation and it is of substantial concern for the eight organic azide compounds in the test set. Indeed, Dorofeeva et al. [21] note that "enthalpies of formation of organic azides are scanty and not always reliable." To address this problem Dorofeeva et al. [21] used a combination of multilevel ab initio model chemistries and density functional theory coupled with the use of isodesmic, isogyric, and other balanced equations to calculate and recommend Δ 298,g values for twenty-nine organic azides including the eight organic azide compounds in our test set. The implication of their work is that their calculated values may be more accurate than the experimental values that were originally available to Byrd and Rice [6,7]. This same group [20,[22][23][24] conducted similar studies of the Δ 298,g values for fifty-seven other energetic nitrogen-containing compounds and recommended the use of several Δ 298,g values that are slightly different than those that were available to Byrd and Rice [6,7].
In light of these findings we compared computed data with the original experimental data set used by Byrd and Rice [6,7] and with a newer alternative experimental reference set (also referred to elsewhere herein as the most recently recommended experimental/reference values available) in which a total of eighteen Δ 298,g reference values recommended by Dorofeeva et al. [18,20,21] were substituted for the original corresponding experimental values in the test set [6,7]. We also note that at least three experimental values      Tables 2 and 3 supports the conclusion that the alternative recommended values, as a group, are indeed more accurate than the original experimental values that were available to and used by Byrd and Rice [6,7] in their investigations. Except for the DFT/Atomic and Group contribution approach the mean absolute error (MAE) was greater for the original experimental reference set than for the new alternative experimental reference set. The fact that the opposite was found for the DFT/Atomic and Group contribution approach is not surprising as many of the compounds in the original experimental reference set were used to parameterize and develop this computational approach. Nevertheless, this approach still performed quite well when compared to the alternative experimental reference set. Our results also suggest that the accuracy of the DFT/Atomic and Group contribution approach might benefit from reparameterization using the alternative experimental reference set.
Because the alternative experimental reference set appeared to be more accurate than the original reference set it was used for the comparisons described below.
The atom and group equivalent DFT approach described by Byrd and Rice [6,7] and the T1 procedure both yielded calculated results in which 51% (23/45) of the predicted Δ 298,g values were within ±2.0 kcal/mol of values in the reference set. The rms deviation from experiment for the predicted T1 gas-phase heat of formation values was found   Dinitromethane + methane →2 nitromethane 2 2 dinitromethane → methane + tetranitromethane 3 Methyl azide + trinitromethane → methane + azidotrinitromethane 4 N-Ethyl-N-nitroethanamine + dimethylamine → N-ethyl-N-ethanamine + N-methyl-N-nitromethanamine (DMNO) 5 1,1-Dinitropropane + methyl azide → ethane + 1-azido-1,1-dinitroethane 6 2 1,1,1-Trinitropropane → butane + hexanitroethane 7 3 propylnitrate →2 propane + nitroglycerin 8 RDX + 3 nitrosyl hydride →3 HNO 2 + TTT (hexahydro-1,3,5-trinitroso-1,3,5-triazine) 9 3 1-nitropiperidine →2 cyclohexane + RDX 10 1-Nitrosopiperidine + TTT →2 1,4-dinitrosopiperazine 11 Piperazine + 2 dimethylnitramine → 2 dimethylamine + 1,4-dinitropiperazine 12 Bis-2,2,2-trinitroethylamine + HNO 2 → H 2 + N-nitro-bis-2,2,2-trinitroethylamine 13 Benzene + 4-nitroaniline → aniline + nitrobenzene 14 Phenol + nitrobenzene → benzene + 2-nitrophenol 15 Phenol + nitrobenzene → benzene + 3-nitrophenol 16 Phenol + nitrobenzene → benzene + 4-nitrophenol 17 Aniline + nitrobenzene → benzene + m-nitroaniline 18 Aniline + nitrobenzene → benzene + p-nitroaniline 19 Methyl azide + benzene → methane + azidobenzene 20 Methyl azide + nitrobenzene → methane + 1-azido-4-nitrobenzene Toluene + nitrobenzene → benzene + PNT (1-methyl-4-nitrobenzene) 23 Toluene + 2 nitrobenzene →2 benzene + 2,4-DNT (1-methyl-2,4-dinitrobenzene) 24 Methyl azide + toluene → methane + azidomethylbenzene 25 Methyl azide + 3-ethylpentane → methane + 3-azido-3-ethylpentane 26 Toluene + 1,3,5-trinitrobenzene → benzene + TNT (2,4,6- to be 2.9 kcal/mol with a mean absolute error (MAE) of 2.3 kcal/mol and a maximum deviation of 6.9 kcal/mol. This compares with an rms deviation of 3.3 kcal/mol, a MAE of 2.4 kcal/mol, and a maximum deviation of 9.9 kcal/mol using data calculated by Byrd and Rice [6,7]. It should be noted that the T1 procedure has been previously shown to reproduce experimental heat of formation values of a large (1805) set of organic molecules with an rms deviation of 2.75 kcal/mol and a MAE of 2.03 kcal/mol [8]. Clearly, the rms and MAE values for the T1 approach reported here agree well with those values reported for the larger data set. These results led us to compare older approaches used to predict/calculate Δ 298,g values. Specifically, we used Benson's group additivity (or group contribution) approach [10,11]. We also used density functional theory coupled with the use of selected isodesmic, isogyric, and other balanced equations. One disadvantage of Benson's group additivity approach is that group values for some of the substructures of several compounds in the test set are not available. We were partially able to overcome this problem by using existing data which allowed us to calculate some of these group values. In total we were able to use Benson's group additivity approach to calculate Δ 298,g values for 36 of 45 compounds in the test set. This approach yielded calculated results in which 64% (23/36) of the calculated Δ 298,g values were within ±2.0 kcal/mol of experimental values. The rms deviation from experiment for the predicted group additivity gas-phase heat of formation values was 5.5 kcal/mol with an MAE of 3.1 kcal/mol with a maximum deviation of 21.2 kcal/mol for HNS. A relatively large deviation from experiment of 14.8 kcal/mol was also found for nitroglycerin. A significant drawback to the use of group additivity theory is that calculated Δ 298,g values are sometimes affected by structural features that are not adequately addressed by theory. Although Δ 298,g values for some strained compounds can be accurately predicted by adding correction factors for structural characteristics such as ring strain, the presence of gauche carbons, and C1-C5 repulsions, contribution of certain other structural features do not appear to be adequately addressed by group additivity theory. HNS and nitroglycerin seem to fall into this category. When these compounds were eliminated from the test set, the rms deviation from experiment for the predicted group additivity Δ 298,g values was found to be 3.5 kcal/mol with an MAE of 2.2 kcal/mol and a maximum deviation of 11.5 kcal/mol for the remaining 34 compounds. The relevant point is that group additivity theory is often quite accurate; however, one must be careful regarding its application.
Density functional theory coupled with the use of selected balanced reactions was also used to calculate Δ 298,g values for comparison to the test set. The reactions used in this investigation are presented in Table 4. We determined that the rms deviation from experiment for the calculated Δ 298,g values was 1.5 kcal/mol with an MAE of 1.1 kcal/mol and a maximum deviation of 3.8 kcal/mol. This approach yielded calculated results in which 82% (37/45) of the Δ 298,g values were within ±2.0 kcal/mol of experimental values. Despite the fact that, for time considerations, we selected a medium size basis set, density functional theory coupled with the use of isodesmic and other balanced equations was still the most accurate approach of the six methods we compared. Although accurate, it is necessary to mention some problems with this approach. First of all, different balanced equations can lead to different Δ 298,g values. Often such differences are small. But sometimes such differences can be substantial. Cramer [12] has noted that "the construction of an isodesmic equation is something of an art, depending on chemical intuition and available experimental data." It should also be mentioned that, when computation time is not a consideration, the use of larger basis sets might be expected to result in even greater accuracy. Indeed, Dorofeeva and colleagues [21][22][23] have have shown that the use of high level computational theory coupled with multiple isodesmic, isogyric, and other balanced equations can result in calculation of Δ 298,g values that are sufficiently accurate that they feel confident in recommending consensus values, in many cases, where there is discrepancy between experimental values and/or between computed values.
Finally, two semiempirical models (RM1 and PM7) were used to calculate Δ 298,g . Semiempirical models are attractive because they are very rapid and their accuracy continues to improve [13,14]. However, RM1 and PM7 were found to be less accurate than the other approaches used to calculate Δ 298,g values for compounds in the test set used for this investigation. Even when Δ 298,g values having very large (greater than 14 kcal/mol) deviations from reference values were omitted, rms and MAE values for calculated RM1 and PM7 data sets were substantially greater than those of the four other computational approaches investigated (Tables 2  and 3).

Conclusion
Because substantial interest exists in developing relatively rapid and accurate procedures for predicting gas-phase heat of formation values for energetic compounds, our overall goal was to determine which of several computational approaches is most suitable for this endeavor. None of the computational approaches were uniformly accurate within ±2.0 kcal/mol of experimental values (i.e., the presumably more accurate alternative experimental values).
However, four of the computational approaches investigated were able to meet this level of accuracy greater than 51% of the time. Moreover, these four computational approaches were accurate to within ±4.0 kcal greater than 77% of the time. For relatively simple compounds in which correction factors can be used to account for ring strain, and so forth, and other factors do not make substantial contribution to Δ 298,g , we conclude that the group additivity approach is about as accurate as those approaches based on higher levels of theory. For compounds that are structurally more complex, the approach using DFT coupled with the use of isodesmic, isogyric, and other balanced equations is more accurate, but more time consuming, than the DFT/Atomic and Group contribution method described by Byrd and Rice [6,7] and the T1 method [8] for this set of compounds. However, since calculated values are often similar and mutually supportive it seems that a reasonable approach would be to use any (or all) of these four procedures when this level of accuracy is acceptable and the goal is to compare predicted properties of new or theoretical compounds with those of structurally similar compounds whose experimental values and/or computed values are well established.