Accurate or Fast Prediction of Solid-State Formation Enthalpies Using

3 hours ago - In view of its success in predicting sublimation enthalpies of molecular crystals near triple point conditions [Ind. Eng. Chem. Res. 201...
0 downloads 0 Views 1MB Size
Subscriber access provided by Kaohsiung Medical University

Thermodynamics, Transport, and Fluid Mechanics

Accurate or Fast Prediction of Solid-State Formation Enthalpies Using Standard Sublimation Enthalpies Derived From Geometrical Fragments Didier André Mathieu Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.8b03001 • Publication Date (Web): 10 Sep 2018 Downloaded from http://pubs.acs.org on September 11, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Accurate or Fast Prediction of Solid-State Formation Enthalpies Using Standard Sublimation Enthalpies Derived From Geometrical Fragments Didier Mathieu∗ CEA, DAM, Le Ripault, 37260 Monts, France E-mail: [email protected] Phone: +33 (0)247344185. Fax: +33 (0)247345158

Abstract In view of its success in predicting sublimation enthalpies of molecular crystals near triple point conditions [Ind. Eng. Chem. Res. 2012, 51, 2814], the geometrical fragment (GF) approach is presently used to estimate standard values ∆sub H 0 adjusted to the reference temperature of 298.15 K. The resulting GF0 scheme is fitted against theoretically confirmed ∆sub H 0 data for 185 high nitrogen compounds and validated using organic compounds not used in the fitting process. It is subsequently used to predict the solid-state formation enthalpy ∆f H 0 (cr) of 40 energetic materials, combined with either correlated ab initio methods or simple models to estimate the gas-phase contribution. The combination of GF0 with ab initio data yields ∆f H 0 (cr) values with state-of-the-art accuracy. Combined with the simple atom pair contribution (APC) model, this procedure predicts ∆f H 0 (cr) with better accuracy and stronger physical grounds than alternative low-cost methods.

1

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1

Introduction

The sublimation enthalpy of a molecular crystal is one of its most fundamental property. It is a valuable indicator of crystal stability, and may be used to evaluate other properties of practical significance such as vapor pressure. 1 Moreover, it is a significant contribution to the solid-state formation enthalpy ∆f H 0 (cr) of organic crystals required to evaluate the energy content, performances and sensitivities of candidate materials of potential interest as components of fuels, propellants or explosives. 2,3 Although ∆f H 0 (cr) may be estimated directly from molecular structure, ie. without explicitly evaluating the gas-phase and sublimation enthalpies, such procedures are deprived of physical bases and rely mostly on empiricism. 4–7 Since the intra- and intermolecular contributions to ∆f H 0 (cr) depend on distinct interactions, they should ideally be evaluated separately. As reviewed recently, 8 a wealth of methods are available for the intramolecular contribution, ie. the gas-phase formation enthalpy. In contrast, a relatively small number of models are available to predict the intermolecular contribution, ie. the sublimation enthalpy. They are commonly classified into four main categories: molecular simulations of the crystal, 9–14 additivity methods, 15,16 methods based on a quantitative structure−property relationship (QSPR) involving so-called General Interaction Properties Function (GIPF) descriptors 17–23 and other QSPR schemes resorting to more general descriptors. 24–28 Numerical simulations are potentially very accurate. However, they have the major disadvantage of requiring knowledge of the crystal structure. For many applications, direct relationships between sublimation enthalpy and molecular structure are needed as crystal data is not available. Many methods reported in current literature to estimate sublimation enthalpies predict ∆sub H(T ) values measured near triple point conditions. 16,25,27,29 Decomposing ∆sub H(T ) into transferable additive contributions associated with so-called geometrical fragments (GFs) proves especially successful, leading to a root mean square error (RMSE) with respect to experiment of only 4 kJ/mol for an external test set obtained from a split of an extended collection of 1300 organic compounds (excluding salts) previously introduced to develop an 2

ACS Paragon Plus Environment

Page 2 of 32

Page 3 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

artificial neural network (ANN) model. 16 In contrast to a subsequent claim, 30 no entries were arbitrarily removed from the dataset on going from the ANN to the GF model. The large errors of about 50 kJ/mol reported for glycine and alanine 31 arise from the incorrect use of the gas-phase structures of these molecules, i.e. respectively H2 N-CH2 -COOH and − H2 N-CH(CH3 )-COOH, instead of the zwitterionic structures H+ and H+ 3 N-CH2 -COO 3 N-

CH(CH3 )-COO− relevant to the crystals. Assuming a value of 90.1 kJ/mol for the contribution of the primary ammonium group −NH+ 3 yields GF estimates of respectively 133.5 kJ/mol and 136.1 kJ/mol, in close agreement with the experimental values of 134.4 kJ/mol and 135.2 kJ/mol. 32 However, it should be stressed that additivity methods are in principle ill-suited for ionic crystals. Indeed, they do not account for the decrease of the magnitude of Coulomb interactions (a major contribution to the cohesion of ionic systems) as the ion sizes increase. 33 A better alternative to the GF model for ionic crystals is provided by volume-based thermodynamics. 34–36 Although satisfactory methods are available to estimate sublimation enthalpies near triple point conditions, the standard value ∆sub H 0 is required in many applications, especially when solid-state formation enthalpies are needed. This requires some adjustments from the mean temperature of measurement to the standard reference temperature T 0 = 298.15 K. In the crudest approximation, ∆sub H 0 can be estimated from ∆sub H(T ) as follows: ∆sub H 0 = ∆sub H(T ) + 2R(T − T 0 )

(1)

where R is the ideal gas constant and T the temperature. This would imply that ∆sub H(T ) decreases little with temperature (by 1.5 K for every 100 K) and differs significantly from ∆sub H 0 only at very high temperatures. Accordingly, the GF approach was recently used to predict heats of compbustion of energetic materials with excellent accuracy. 37 Unfortunately, the large differences actually observed for some compounds between sublimation enthalpies measured at various temperatures clearly show that this approximation

3

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

may underestimate the temperature dependence of sublimation enthalpies. 38 This supports the view that the 2R approximation sometimes "introduces perhaps more uncertainty than accuracy", as pointed out by Maschio et al. 13 Therefore, the latter cannot be taken for granted nor considered as a safe approach to derive ∆sub H 0 from accurate ∆sub H(T ) data. In fact, despite the success of the GF approach for sublimation enthalpies at triple point conditions, a satisfactory approach for the standard value ∆sub H 0 is still lacking. Among current models, only two have been extensively validated for general organic compounds. The first one is a QSPR based on so-called HYBOT descriptors, 26 and the second one is a group contribution (GC) method. 30 The former requires a proprietary software package whose details are not fully disclosed in literature, whereas the latter exhibits other issues, as detailed below. In this context, the best demonstrated performance is currently obtained within the above-mentioned GIPF approach. 22,31 The GIPF descriptors are derived from the molecular electrostatic potential (MEP) computed on an isoelectronic density surface surrounding the molecule. They are defined so as to reflect major intermolecular interactions like hydrogen bonds and Coulomb interactions. The beauty of this approach introduced many years ago by Murray et al. 19 stems from the fact that it requires only 3 fitting parameters, making it possible to derive a GIPF model from a small dataset. Using a database of reliable experimental values confirmed by high-level theoretical calculations for 185 high-nitrogen compounds (hereafter denoted SUB-185), Suntsova and Dorofeeva (SD) recently introduced a modified GIPF model (mGIPF) involving an additional fitting parameter, which stands as the current state of the art for standard sublimation enthalpies. 31 Nevertheless, GIPF-based models have a number of disadvantages compared to the GF approach, including a much higher computational cost and a spurious quadratic dependence of the sublimation enthalpy on the molecular surface area that prevents its application to compounds very different in sizes. 18 Due to their computational requirements, GIPF models have not yet been assessed against extensive datasets. In this context, a well-validated GFlike model for ∆sub H 0 is desirable. The present article reports such a method. To emphasize 4

ACS Paragon Plus Environment

Page 4 of 32

Page 5 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

the fact it is fitted against standard values ∆sub H 0 rather than values measured at near triple point conditions, it is hereafter referred to as the GF0 scheme.

2 2.1

Modeling approach and datasets Training set

Compared to the previously reported GF procedure for ∆sub H(T ) data measured at triple point conditions, 16 GF0 mainly differs in the fact that it is fitted against standard values of the sublimation enthalpies. The fact that both the original and modified GIPF equations considered by Suntsova et al. (SD) 31 perform remarkably well suggests that this good performance owes to the high quality of the SUB-185 dataset. The latter is therefore used in the present work to fit GF0 . As a result, there are two reasons why GF0 should be better suited than the earlier GF method 16 to design energetic materials, as SUB-185 is 1) made only of high-quality data, and 2) focused on high nitrogen compounds.

2.2

Model equation

Because of the reduced size and diversity of SUB-185 compared to the training set used in ref 16, the number of fitting parameters is decreased from 35 to 18 on going from GF to GF0 . Furthermore, the model is restricted to CHNO compounds and earlier crowding corrections are irrelevant. 16 Finally, in addition to fragment contributions hk , only three ring corrections are presently used, namely Ra for any aromatic ring, Rs for small aliphatic rings (with only 3, 4 or 5 atoms) and Rl for any larger ring. Accordingly, ∆sub H 0 is obtained as follows: ∆sub H 0 =

X

nk hk + na Ra + nl Rl + ns Rs

(2)

k

where index k runs over all geometrical fragments, ie. any non-hydrogen atom with corresponding hydrogen neighbors, nk is the number of fragments of type k in the molecule, 5

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

na , nl and ns are the number of rings associated with Ra , Rl and Rs , respectively. The values of these parameters were derived from an ordinary least squares regression using the statsmodels package. 39 For clarity, worked out examples of application of GF0 are provided in section S1 (Supporting Information). Since the present number of fitting parameters remains large compared to GIPF methods, GF0 is fitted against the whole SUB-185 dataset. It is then validated through the three following procedures. First, a leave-one-out (LOO) cross-validation (CV) against the training set is carried out. Secondly, the model is applied to external test sets denoted SUB-095 and SUB-996. Thirdly, ∆sub H 0 values estimated using the present GF0 scheme are further validated for energetic materials through their use in evaluating solid-state formation enthalpies for an additional dataset (SOL-40) already used to test the mGIPF model. 31

2.3

Applicability domain

Just like any method based on a simple analytical equation, the GF model is likely to lead to large errors in some specific cases. Obvious limitations of the approach include: 1. The complete neglect of the electron distribution. In particular, the contribution of any atom to ∆sub H 0 does not depend on the surrounding electron density; 2. The fact that any atom contribution is assumed to depend only on local features, ignoring the overall structure of the molecule; 3. The fact that molecular geometry and 3D conformation are ignored. As a consequence, non-bonded intramolecular interactions are neglected, although they should decrease ∆sub H 0 as the availability of the interacting groups for intermolecular interactions is decreased; 4. The assumption that the accessibility of every atom depend only on its first neighbors, ignoring the role of the molecular shape;

6

ACS Paragon Plus Environment

Page 6 of 32

Page 7 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

5. The purely additive character of the method, which is ill-suited to describe the subtle cooperativity effects driving the structure and energetics of crystals with extensive hydrogen bonded networks. In view of these approximations, some molecules clearly lie outside the scope of the GF model. They fall into two main categories. The first one includes compounds which cannot form a close-packed crystal owing to their rigid shape (like fullerenes and their derivatives) or to the presence of voids or concavities that cannot be filled by neighboring molecules, as exemplified by [1.8]paracyclophane (Fig. 1). Due to the central void, carbon atoms in this compound interact with surrounding molecules only on one side, while the model assumes a close-packed structure. Consequently, the current GF0 model would overestimate the measured ∆sub H 0 value of 110.9 kJ/mol by as much as 74 kJ/mol for this specific compound.

Figure 1: Van der Waals model of [1.8]paracyclophane (C21 H26 ): 3D structure obtained by global optimization of a starting random conformer described using the MMFF94 force field through 1000 iterations of the simulated annealing algorithm implemented in the TINKER package. 40 In the lack of 3D molecular models, it is difficult to reliably predict the occurrences of such crystal cavities. In this work, it is assumed that they are likely to be encountered for large molecules (over 20 carbon atoms) with at least one set of three adjacent rings (i.e. 7

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 32

sharing at least one common atom) including at least one aliphatic (to avoid including fused aromatic systems like polyaromatic hydrocarbons in this definition). Such molecules are therefore excluded from the applicability domain of GF0 . The second category of excluded molecules are carbohydrates, as these compounds are well-known to potentially exhibit complex hydrogen bond networks in crystals. 41 In this article, carbohydrates are defined as saturated molecules with empirical formula Cm H2n+p On where p ≤ 2. Although p = 0 for most carbohydrates, especially monosaccharides, positive values of this parameter allow to encompass less-common sugars like deoxyribose.

2.4

External test sets

As mentioned above, two external test sets of ∆sub H 0 data are used for a direct validation of GF0 . The first one (SUB-095) includes 95 molecules from the set of 158 compounds studied by McDonagh et al. 42 The second one (SUB-996) include 996 compounds from the database compiled by Acree and Chickos, hereafter referred to as the AC database. 38,43 These two test sets are obtained from the original datasets after compounds beyond the applicability domain of GF0 or already present in SUB-185 have been removed.

2.5

Validation through solid-state formation enthalpies

Since only experimental values of ∆f H 0 (cr) are available for the 40 high nitrogen compounds in SOL-40, these data only allow an indirect validation of present ∆sub H 0 predictions through the equation: 0 0 ∆f Hsol = ∆f Hgas − ∆sub H 0

(3)

0 where theoretical values are used for the gas-phase formation enthalpy ∆f Hgas . The latter

are either high-level ab initio values computed and the G4 44 or G4MP2 45 level of theory and reported in ref 31, semi-empirical PM3 values 46 presently computed using the Firefly software package, 47 or even simpler estimates derived from the recent Atom Pair Contribution (APC) 8

ACS Paragon Plus Environment

Page 9 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

additivity scheme. 8

3

Model parameters

The fit of GF0 against SUB-185 yields 18 parameters whose values are listed in Table 1. The associated standard deviations δi reveal that two parameters, namely C2-0 and N1-0, are statistically ill-defined. This is because the two corresponding fragments are usually involved in −CN groups. As a result, their number is the same for every molecule presently considered, except for triphenylmethyl azide where the terminal =N atom belongs to a −N3 group. It is possible to get all parameters well-defined by dumping the C2-0 and N1-0 parameters into a single contribution associated with the cyano group. 48 However, this would make the present model unapplicable to azides. In fact, the values presently obtained for these two parameters are quite reasonable. As shown on Fig. 2, the same parameters would be obtained if N1-0 was linearly extrapolated from N2-0 and N3-0. Therefore, although the individual values of these two parameters effectively depend on the data for a single molecule, it may be assumed that their combined use does not entail any significant error. Fig. 2 and Fig. 3 emphasize the fact that the contribution of a given atom to ∆sub H 0 varies smoothly as its number of hydrogen neighbors or total number of neighbors are varied. More specifically, Fig. 2 illustrates the steady decrease of the contribution hk associated with fragment k as the coordination number of the central atom increases, due to the fact that intermolecular interactions get increasingly hindered by the surrounding neighbors. Remarkably, although hk depends significantly on the total coordination number of the atom, its value exhibits only a weak dependence on the atomic number as long as only elements C-N-O are considered. Fig. 3 shows how hk increases on substituting neighbors of the central atom with hydrogen atoms, as a result of two factors: the increased availability of this atom arising from the fact that a neighbor is substituted with a much smaller atom, and the fact that the contribution

9

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 1: Additive contributions hi to standard formation enthalpies ∆sub H 0 with associated standard deviations δi and numbers of occurences Niocc . For clarity, atomic fragments are shown with bonds and attached hydrogen atoms made explicit, as done in a recent paper 49 i C4-3 C4-2 C4-1 C4-0 C3-1 C3-0 C2-0 N3-2 N3-1 N3-0 N2-0 N1-0 O2-1 O2-0 O1-0 Ra Rl Rs

atom −CH3 >CH2 >CH− >C< =CH− =C< ≡C− −NH2 >NH >N− =N− ≡N −OH >O =O

hi δi Niocc 24.0 1.8 50 5.5 1.3 36 -27.7 4.1 5 -48.2 5.1 13 3.6 0.7 137 -16.0 2.3 155 2.5 12.6 11 40.3 1.8 36 19.7 2.1 46 -16.3 3.1 92 6.8 1.3 65 29.1 12.8 12 37.2 2.4 24 4.1 2.9 12 26.0 2.1 124 40.8 5.1 141 29.7 4.2 13 43.2 4.5 17

of hydrogen atoms are dumped into that of the central atom. As a consequence, hk increases by about 20−25 kJ/mol. The only exception (red arrow on Fig. 3) is the large increase by 36 kJ/mol observed for nitrogen on going from >N− to >NH (i.e. from N3-0 tu N3-1) as a consequence of hydrogen bonding. Similar trends were observed in the previous application of the GF approach to ∆sub H(T ) 16 and for other properties such as densities 50,51 or solubilities. 49 This possibility to check the physical consistency of the model parameters is a significant advantage of GF models over more heavily parametrized GC methods. Together with their relatively small number, the interpretability of the GF parameters enables the derivation of predictive models on the basis of smaller datasets than possible using standard GC methods.

10

ACS Paragon Plus Environment

Page 10 of 32

Cx-0 Nx-0 Ox-0

20

0

-20

-40

1

2 3 Total coordination number

4

Figure 2: Fragment contributions to ∆sub H 0 plotted according to the total coordination number of the central atom, for fragments deprived of hydrogen atoms.

40 Contribution to sublimation enthalpy (kJ/mol)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Contribution to sublimation enthalpy (kJ/mol)

Page 11 of 32

20

0

-20 C3-x C4-x N3-x

-40

0

1 2 Number of hydrogen neighbors

3

Figure 3: Fragment contributions to ∆sub H 0 plotted according to the number of hydrogens attached to the central atom.

11

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

4

Results

All ∆sub H 0 values estimated using the present model are provided as SI together with corresponding experimental values. These data are reported in Table S1 for SUB-185, in Table S2 for SUB-095 and SUB-996, and in Table S3 for SOL-40. Predicted values of ∆sub H 0 are plotted against experiment on Fig. 4. For any compound in SUB-996, the horizontal error bar indicates the ranges of values spanned by the experimental data in the AC compilation, 38,43 while the symbol represents the most recent value, which is also the value used for the calculation of statistical data reported in Table 2 where the performance of GF0 for the various datasets is summarized. SUB-185 (training set) SUB-095 SUB-992

250

0

∆subH calculated (kJ/mol)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 32

200 U SV

150 B C A

R

100 D

W

50 Y

Z T

0

0

50

X

100 150 200 0 ∆subH observed (kJ/mol)

250

Figure 4: Sublimation enthalpies: present GF estimates versus AC data. Major outliers are as follows: A: hexanitroethane C2 N6 O12 ; B: [2.2]perhydroparacyclophane C16 H28 ; C: [2.2]metaparacyclophane C16 H16 ; D: nortricyclene C7 H10 ; Z: (5α)-cholestane C27 H48 ; Y: 1,1’Diadamantyl C20 H30 ; X: Cookson diketone C11 H10 O2 ; W: 1,6-anhydro-β-D-glucose; V: 1docosanol C22 H46 O; U: tetracosanoic acid C24 H48 O2 ; T: isogarudane C14 H16 ; S: dibenzo-24crown-8 ether C24 H32 O8 ; R: pentacene C22 H14 . 12

ACS Paragon Plus Environment

Page 13 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

4.1

Training set

The RMSE for SUB-185 is only 10.4 kJ/mol, as expected from the significant number of fitting parameters involved. Calculated data lie mostly within 25 kJ/mol from experiment, the only notable exception being for 3-hydroxypyridine N-oxide for which the calculated enthalpy is lower than experiment by about 36 kJ/mol. The most serious underestimations are generally for N-oxides, with deviations from experiment as large as −26 kJ/mol for isonicotinamide N-oxide, −21 kJ/mol for quinazoline 1,4-dioxide or −20 kJ/mol for 4nitrobenzylidene tert-butylamine N-oxide. This is understandable as the specially negative charge on the corresponding oxygen atoms is ignored by the present model. Other data for N-oxides are in very good agreement with experiment, or even overestimated by +11 kJ/mol for 2,4,6-trinitropyridine N-oxide. The latter case can be rationalized by invoking the electron withdrawing effect of the three nitro groups on the aromatic ring which is expected to decrease the polarity of the dative bond. On the other hand, the most overestimated values are for 1,4-diazabicyclo[2.2.2]octane (+28 kJ/mol) and styphnic acid (+26 kJ/mol). It should be kept in mind that in contrast to the usual practice in QSPR approaches, the leave-one-out (LOO) cross-validation (CV) is not used here to to tune some model hyperparameters nor to select an optimal model. Therefore, the LOO-CV statistics reflects not only the robustness of the model, but also its predictive value for high nitrogen compounds similar to those included in the training set. With a corresponding RMSE of 11.9 kJ/mol, the performance of the present model is better than for the original GIPF model (RMSE ' 15 kJ/mol) 52,53 and similar to that obtained using the recent GIPF and mGIPF models (RMSE=11.2-11.7 kJ/mol) developed by SD on the basis of the same data set. 31 In other words, despite its simplicity, GF0 provides state of the art accuracy with regard to evaluating ∆sub H 0 of high nitrogen compounds without resorting to crystal data.

13

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

4.2

Page 14 of 32

Test sets

Not surprisingly, GF0 predictions are less accurate for general organic compounds as they may differ significantly from the high nitrogen compounds in the training set. Good results are still obtained as long as only simple moderate-size molecules are considered, as clear from the RMSE of only 13.5 kJ/mol obtained for SUB-095, and the fact that all predictions are within 35 kJ/mol from experiment for this dataset. This performance is close to that obtained using advanced QSPR models specifically developed against this dataset, with RMSE values ranging from 10.7 kJ/mol for partial least square regression to 11.9 kJ/mol for support vector regression. 42 Table 2: Fitting ability and predictive value of GF0 as characterized by the values of the determination coefficient (R2 ), average absolute deviation (AAD) and root mean square error (RMSE). The minimal and maximal deviations from experiment are reported as MIN and MAX, respectively. AAD, RMSE, MIN and MAX in kJ/mol. Dataset SUB-185 SUB-185 SUB-095 SUB-996

status training LOO-CV test test

N 185 185 95 996

R2 0.64 0.53 0.63 0.69

AAD 8.0 9.2 10.8 13.5

RMSE 10.4 11.9 13.5 18.4

MIN -36 -38 -34 -100

MAX +28 +31 +24 +48

Application of GF0 to the extended test set SUB-996 reveals more significant errors, including 17 molecules with errors > 50 kJ/mol, in all cases because GF0 underestimates ∆sub H 0 . Consequently, a significantly larger RMSE value of 18.4 kJ/mol is obtained for SUB-996. Despite the fact that the AC database is carefully curated and up-to-date, significant experimental errors cannot be ruled out. For instance, this database includes a value of 226.8 kJ/mol for the standard sublimation enthalpy of the CL-20 explosive (CASRN: 135285-90-4), whereas a value of only 157.6 kJ/mol is recommended by SD. 31 This demonstrates that deviations from experiment as large as 70-100 kJ/mol are to be expected when trying to predict ∆sub H 0 from molecular structure. For most compounds with significant differences between calculated and measured values, a single experimental value is available. 14

ACS Paragon Plus Environment

Page 15 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

This is the case for (5α)-cholestane, for which the most serious error is observed. Although the only experimental value of 133.8 kJ/mol available for this molecule has not been subsequently confirmed, the GF0 value of 34.1 kJ/mol calculated for this large molecule (C27 H48 ) is unrealistically small. This is a consequence of the large number of >CH− fragments in this molecule, for which the associated negative contribution to ∆sub H 0 over-corrects the fact that for each non-hydrogen neighbor of the central carbon, intermolecular interactions are hindered by the presence of at least two geminal non-hydrogen atoms. For 3,4-dihydroxy3-cyclobutene-1,2-dione, two widely conflicting values are reported. The calculated value of 103.5 kJ/mol is much lower than the retained experimental value of 154.3 kJ/mol reported in 1983, 54 but fairly consistent with the value of 83.7 ± 16.7 kJ/mol derived from earlier measurements. 55,56 In spite of experimental uncertainties, many errors apparent on Fig. 4 clearly reflect limitations of GF0 , and presumably of any current predictive scheme based only on 2D structural formulas. Notwithstanding (5α)-cholestane, severely underestimated enthalpies are mostly for compounds with small cage structures (1,1’-diadamantyl, Cookson diketone, isogarudane...) It might be tentatively assumed that the specially high number density of atoms in such structures makes it possible for any atom to interact with a specially large number of other atoms from neighboring molecules within a relatively short distance, hence contributing to enhanced sublimation enthalpies. Other significant deviations from experiment are easier to explain. The fact that ∆sub H 0 is underestimated for 1,6-anhydro-β-D-glucose is clearly related to the ability of this sugar to form extended hydrogen bond networks with cooperative effects, just like in carbohydrates which are excluded from the GF0 applicability domain. The largest overestimation (+48 kJ/mol) is for hexanitroethane, a crowded molecule in which the interactions of a nitro group with the molecular surroundings is specially weak due to steric hindrance. A significant overestimation (+35 kJ/mol) is also observed for trinitromethane, presumably for similar reasons. 15

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Notwithstanding hexanitroethane and trinitromethane, significantly overestimated values of ∆sub H 0 are mostly for small cyclophanes, due to the fact that non-bonded atoms from the same molecule are constrained to lie close to each other, thus decreasing the magnitude of intermolecular interactions that contribute to ∆sub H 0 . Larger errors would be obtained for some larger cyclophanes like the one shown on Fig. 1 if such molecules were not excluded from the applicability domain of the model. Although such errors may be intuitively anticipated, it is not obvious to develop a systematic method to identify such outliers automatically. It is worth noting that many of the most serious deviations from experiment are for hydrocarbon compounds, like cyclophanes or cholestane, due to the role of molecular size and shape. RMSE values calculated for the 133 hydrocarbon and 863 non-hydrocarbon compounds in the present dataset are respectively 20.4 and 18.0 kJ/mol, demonstrating that ∆sub H 0 is by no means easier to predict for the former category. Considering average relative errors, they are even larger for hydrocarbons (16.6%) than for other compounds (11.7%). More generally, the reliability of ∆sub H 0 predictions does not significantly depend on the functional groups present in the molecules studied, as observed previously using a GC method. 30 Furthermore, since carbohydrates have been a priori discarded, the same typical accuracy (RMSE=18.5 kJ/mol) is obtained for the 547 molecules with labile protons and for the 449 without. The most obvious determining factor for the reliability of GF0 is molecular size. Sublimation enthalpies are poorly predicted for the 67 molecules with at least 20 carbon atoms (RMSE=28%). However, average relative errors are close to 12% whatever the size of the molecules considered. Finally, it is interesting to compare GF0 with earlier general predictive models for standard sublimation enthalpies. The RMSE value derived from HYBOT-based predictions for a test set of 450 molecules is only 14.9 kJ/mol, instead of 18.4 kJ/mol for GF0 applied to SUB-996. 26 However, this validation was carried out for relatively simple compounds. For instance, the test set considered to assess this earlier model included only one compound among the 17 substances from SUB-996 with GF0 errors > 50 kJ/mol. Similarly, the low 16

ACS Paragon Plus Environment

Page 16 of 32

Page 17 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

RMSE of 10.8 kJ/mol obtained for a dataset of 1269 organic compounds using the GC model of Gharagheizi et al. is remarkable. 30 However, despite the large training set used (1015 entries), this model is specially prone to overfitting issues. Indeed, as many as 147 group contributions were introduced, about a quarter of which (34) were not represented in the test set. In other words, these parameters make it possible for all training set data to fit into the model, but their ability to positively contribute to predictions for an external test set remains unknown. Furthermore, many substructures are under-represented in the data set. Twelve substructures appear only once in the training set reported by the authors, and three of them are never encountered, which indicates that additional entries were actually used in the fit. Although GF0 may also involve statistically ill-defined parameters, the systematic variations observed on modifying atomic environments (shown on Fig. 2 and Fig. 3) give confidence in the values of the parameters. Furthermore, more heavily parametrized fragment-based methods are unlikely to correct for the most serious limitations presently observed for GF0 . Indeed, problems encountered with cyclophanes or carbohydrates are associated with the ability of 3D conformers to form close-packed crystals or with cooperativity effects determined by the detailed structures of hydrogen bond networks. Such effects are inherently non-additive, as the contribution of any substructure in the molecule significantly depends on the other substructures with which it can interact.

5

Prediction of solid-state formation enthalpies

In this section, the reliability of GF0 for energetic materials is estimated on the basis of its ability to predict the solid-state formation enthalpy of molecular crystals in SOL-40, a dataset previously used to assess a combination of the mGIPF method with high level G4/G4MP2 calculations. 31 In what follows, the 40 crystals considered are referred to using the same indexes or identifiers as used in Table 3 of ref 31. This data is duplicated in Table S3 for convenience.

17

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

In principle, compound 6 (1-methyl-2-nitroguanidine) lies beyond the scope of GF0 due to a lack of the N2-1 contribution to be associated with the =NH fragment on this molecule. Because this fragment exhibits a labile proton bounded to a nitrogen atom, we arbitrarily assigned to this fragment a contribution of 40 kJ/mol similar to that obtained for −NH2 . This value might be expected to be slightly underestimated as it does not account for the enhanced availability of the lone pair on the N atom.

5.1

Results obtained using ab initio gas phase data 100

0 Calc. - exp. (kJ/mol)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 32

-100

-200

-300 Suntsova & Dorofeeva 0 G4/G4MP2 - GF

-400 0

5

10

15 20 25 Molecule index

30

35

40

Figure 5: Deviation of calculated solid state formation enthalpies from experiment for 40 high nitrogen compounds. Grey circles represent the results of recent state-of-the-art calculations reported in ref 31, and black squares are for present results derived from GF0 values of the sublimation enthalpy and G4 or G4MP2 data reported in 31 for the gas-phase contribution. Theoretical values of ∆f H 0 (cr) calculated on the basis of ab initio calculations at the G4 or G4MP2 levels are first considered. Fig. 5 shows how deviations from experiment previously reported in ref 31 are affected as mGIPF is substituted with GF0 for the determination of ∆sub H 0 . In abscissa is reported the molecule index as originally defined in ref 31. There are only three compounds with deviations > 50 kJ/mol, namely 16, 30 and 32. The enthalpies calculated for 16 and 30 are very close to previously reported mGIPF data, and much lower than experimental values which are presumably incorrect. 31 For 32, the 18

ACS Paragon Plus Environment

Page 19 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

GF0 value is lower than experiment by 59 k/mol, but only 23 kJ/mol lower than the mGIPF value. These results reflect the close proximity between mGIPF and GF0 results. According to the previous analysis, 31 the experimental values are especially dubious for 16, 17 and 30. Ignoring these three compounds, the average absolute deviation (AAD) between calculations and experiment is 19.5 kJ/mol for both procedures. This further supports the above mentioned finding that GF0 and mGIPF are of comparable accuracy. This result is remarkable as GF0 is a virtually no-cost procedure, in contrast to mGIPF which requires extensive numerical computations to determine the electronic structure and the corresponding average descriptors based on an isolectronic surface. However, this is no practical advantage if these procedures to estimate ∆sub H 0 are associated with high level ab initio theories like G4 or G4MP2 for the gas-phase contribution, as the determination of 0 is by far the most demanding step in this context. ∆f Hgas

5.2

Results obtained using PM3 or APC gas phase data

It does not harm to resort to GF0 to go from ab initio gas-phase enthalpies to solid-state values, possibly along with mGIPF calculations in order to check the consistency of these two models for ∆sub H 0 . However, in view of its simplicity, GF0 should be especially valuable in the context of the fast screening of candidate molecules for specific applications, especially in the field of high energy materials. 57 Therefore, it is interesting to investigate the accuracy that may be obtained by combining GF0 with faster procedures like APC or PM3 for the gas-phase contribution. The deviations from experiment of solid-state formation enthalpies obtained by substracting GF0 sublimation enthalpies from either APC or PM3 are shown on Fig. 6, where the notation is the same as on Fig. 5. Both models provide fair results for compounds 1-15. The large differences between measured and estimated values of respectively +104.7, −50.1 and +447.7 kJ/mol observed for compounds 16, 17 and 30 and the APC-GF0 level are fairly similar to the corresponding deviations of +72.7, −90.9 and +425.2 kJ/mol reported in ref 31, which confirms the doubts 19

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

400

200 Calc. - exp. (kJ/mol)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 32

0

-200

0

-400

-600

APC-GF 0 PM3-GF

0

5

10

15 20 25 Molecule index

30

35

40

Figure 6: Deviation of calculated solid state formation enthalpies from experiment for 40 high nitrogen compounds. Calculations use GF0 for the sublimation enthalpy, and either APC or PM3 for the gas-phase contribution. regarding experimental data for these three materials. The other materials included in this test set are mostly tetrazines (compounds 21-28 and 31) and furazans (31-40). While experimental measurements are often inaccurate for such unstable compounds, ∆f H 0 (cr) data estimated on the basis of the APC or PM3 clearly exhibit large deviations from experiment (Fig. 6) compared to the values derived from G4 calculations (Fig. 5, note the different scale compared to Fig. 6). At the APC level, the largest differences between calculated and observed data are +249 kJ/mol for 31, +232 kJ/mol for 38 and +155 kJ/mol for 28. It may be observed that 28 and 31 are the two only compounds with the 1,2,3,4-tetrazine-1,3dioxide moiety. This confirms the need to incorporate more specific ring corrections into the APC procedure. 8 After removing 16, 17 and 30, the AAD increases from 19.5 kJ/mol to 59 kJ/mol on going from G4 to APC for the gas-phase contribution, and to 84 kJ/mol if PM3 is used instead. PM3 overestimates ∆f H 0 (cr) by 346 kJ/mol for 40 and by 261 kJ/mol for 32. This confirms the superiority of APC over standard semi-empirical quantum chemical methods observed for general organic compounds as well as more specific families. 8

20

ACS Paragon Plus Environment

Page 21 of 32

5.3

Zero-cost APC-GF0 procedure

To get further insight into the performance of the zero-cost procedure combining APC gasphase enthalpies and GF0 sublimation enthalpies (APC-GF0 ), calculated values and corresponding deviations from experiment are plotted against measured values on Fig. 7. By far the most significant deviation from experiment is obtained for 30, a compound with two s-triazine rings bisubstituted by two azide groups and for which the calculated ∆f H 0 (cr) value is almost 450 kJ/mol lower than the overestimated experimental value of 2171 kJ/mol. As emphasized by the blue arrow, this very large deviation is fixed as the experimental value is substituted with the high level theoretical value reported in ref 31. In view of the large experimental inaccuracies observed for this dataset and especially for 30, inherent to technical issues with the characterization of unstable compounds, quickly evaluating ∆f H 0 (cr) using the APC−GF0 model provides values even better than experiment, although clearly not as reliable as those obtained on the basis of ab initio thermochemical recipes.

Calculated data / deviation from experiment

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Comparison against experiment Comparison against SD data Deviation from experiment

2000

1500

1000

500

0

-500 -500

0

500 1000 1500 reference data

2000

Figure 7: Calculated APC-GF0 solid state formation enthalpies for 40 high nitrogen compounds plotted against reference values (circles) and corresponding deviations (squares). Reference values are either experimental (light green symbols) or high level theoretical data from ref 31, referred to as SD data (dark blue symbols). Notwithstanding 30, other significant deviations are −232 kJ/mol for 38, a large com21

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

pounds with four furazan rings linked together through azo linkages, and −250/−155 kJ/mol for 31/28, two compounds with the 1,2,3,4-tetrazine-1,3-dioxide moiety which is of special interest as a building block for green energetic materials with high performance. 58 Furthermore, 33 is very similar to 31 except for the oxidation of the furazan moiety into a furoxan group. While no experimental data is available for this compound, it is remarkable that the APC value is in error by −252 kJ/mol with respect to the G4/G4MP2-mGIPF value reported in ref 31, just like for 33. This supports the assumption that introducing additive correctiond for ring strains should significantly improve the reliability of the APC model. 8 Further evidence for this is provided by the fact that other significant errors (> 50 kJ/mol) on ∆f H 0 (cr) are systematically observed for compounds with unusual ring structures (furazan and 1,2,4,5-tetrazine), whereas the largest error for an acyclic compound is the overestimation by 50 kJ/mol obtained for 1-methyl-2-nitroguanidine (6).

5.4

Value of GF0 for the design of high energy materials

In view of applications in the field of energetic materials, it is desirable that the solid state formation enthalpy reported on a per mass basis be predicted to within 0.5 kJ/g of experiment. 20 For the present high nitrogen compounds, and ignoring 16, 17 and 30 for which experimental data are probably erroneous, all predictions are within this error range if GF0 is combined with G4/G4MP2 calculations. On going from ab initio to APC data for the gas-phase formation enthalpies, larger deviations are observed for only five compounds, for which this quantity is too low by 0.6−1.6 kJ/g. Therefore, APC-GF0 may be useful to discard candidate compounds that are clearly unsuitable for applications as energetic materials, provided that a safety margin of about 2 kJ/g is taken into account. Since the enthalpy per mass unit ranges between -2 kJ/g to over +6 kJ/g for present materials, this procedure would already be quite effective as a preliminary rapid screening tool, although a better accuracy is necessary for final selection of the most promising candidates, as obtained using ab initio thermochemical recipes instead of APC. A simple script implementing the present 22

ACS Paragon Plus Environment

Page 22 of 32

Page 23 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

GF0 model is provided as Supporting Information to this paper. In view of fast evaluation of solid state formation enthalpies, this script can be used in conjunction with the former one implementing the APC model provided elsewhere. 8

6

Conclusion and perspectives

The present study confirms preliminary findings that ∆sub H 0 is more difficult to estimate than genuine experimental values ∆sub H(T ) measured at near triple point conditions. In particular, while geometrical fragments proves very successful at predicting the latter values, significantly larger uncertainties are noted when applying this approach to ∆sub H 0 . They appear to be due to limitations inherent to additivity methods, especially failures of some compounds to form close-packed crystals and cooperative effects in crystals with extensive hydrogen bonding, like carbohydrates. Nevertheless, GF0 (the present application of geometrical fragments to ∆sub H 0 ) proves competitive with much more costly procedures like those based on the molecular electrostatic potential. In view of its good accuracy and zero cost, this model is recommended whenever standard sublimation enthalpies are needed in the lack of structural data for the crystals studied. If more costly methods can be afforded, GF0 estimates are nevertheless useful as a consistency check. Combined with fast procedures (including the APC model) to estimate gas-phase formation enthalpies and densities, the method should be especially useful for preliminary screening of potential candidate compounds to be used as fuels, explosives or propellants. In a subsequent step, refined estimates of the solid-state formation enthalpy may be obtained for the most promising molecules using composite ab initio methods or other high level quantum chemical calculations for the gas-phase contribution. At this stage, GF0 may still prove useful, although its computational advantage over methods based on the molecular electrostatic potential becomes irrelevant in this context and numerical procedures based on predicted

23

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

crystal structures might be preferred for better accuracy. Further development of the present method, especially beyond CHNO crystals, clearly depends on the availability of high quality databases similar to the one recently reported by Suntsova and Dorofeeva for high nitrogen compounds. In addition, progress regarding quantum chemical calculations in molecular crystals will make it possible to compile accurate lattice energies providing additional data to fit improved analytic relationships. 12–14 Finally, it must be kept in mind that this approach, just like any alternative scheme based on simple analytic relationships, is deemed to fail for open crystals in which molecules are not closely packed, as emphasized in the introduction to this article.

Acknowledgement This work arose from fruitful discussions with Julien Glorian from Institut Saint-Louis, France.

Supporting Information Available The following files are available free of charge. • PDF: Section S1 providing worked-out examples to illustrate the GF0 model • XLSX: Tables S1-S3 with full details of present results • TXT: Python script based on the RDKit library implementing the present model

References (1) Ruz, V.; González, M. M.; Winant, D.; Rodríguez, Z.; Van den Mooter, G. Characterization of the Sublimation and Vapor Pressure of 2-(2-Nitrovinyl) Furan (G-0) Using

24

ACS Paragon Plus Environment

Page 24 of 32

Page 25 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Thermogravimetric Analysis: Effects of Complexation with Cyclodextrins. Molecules 2015, 20, 15175–15191. (2) Mathieu, D.; Beaucamp, S. Matériaux Energétiques; Encyclopédie "Techniques de l’Ingénieur", Traité "Sciences Fondamentales", Fascicule AF 6710: Paris, 2004. (3) Rice, B. M.; Byrd, E. F. C.; Mattson, W. D. High energy density materials; Structure and Bonding; Springer-Verlag: Berlin, 2007; Vol. 125; pp 153–194. (4) Salmon, A.; Dalmazzone, D. Prediction of Enthalpy of Formation in the Solid State (at 298.15 K) Using Second-Order Group Contributions - Part 1. Carbon-hydrogen and Carbon-Hydrogen-Oxygen Compounds. J. Phys. Chem. Ref. Data 2006, 35, 1443– 1457. (5) Salmon, A.; Dalmazzone, D. Prediction of Enthalpy of Formation in the Solid State (at 298.15 K) Using Second-Order Group Contributions - Part 2: Carbon-Hydrogen, Carbon-Hydrogen-Oxygen, and Carbon-Hydrogen-Nitrogen-Oxygen Compounds. J. Phys. Chem. Ref. Data 2007, 36, 19–58. (6) Guella, S.; Argoub, K.; Benkouider, A. M.; Yahiaoui, A.; Kessas, R.; Bagui, F. Artificial Neural Network-Group Contribution Method for Predicting Standard Enthalpy of Formation in the Solid State: C–H, C–H–O, C–H–N, and C–H–N–O Compounds. Int. J. Thermophys. 2015, 36, 2820–2832. (7) Jafari, M.; Keshavarz, M. H.; Noorbala, M. R.; Kamalvand, M. A Reliable Method for Prediction of the Condensed Phase Enthalpy of Formation of High Nitrogen Content Materials through their Gas Phase Information. ChemistrySelect 2016, 1, 5286–5296. (8) Mathieu, D. Atom Pair Contribution Method: Fast and General Procedure To Predict Molecular Formation Enthalpies. J. Chem. Inf. Model. 2018, 58, 12–26.

25

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(9) Arnautova, E.; Zakharova, M.; Pivina, T.; Smolenskii, E.; Sukhachev, D.; Shcherbukhin, V. Methods for Calculating the Enthalpies of Sublimation of Organic Molecular Crystals. Russ. Chem. Bull. 1996, 45, 2723–2732. (10) Beaucamp, S.; Bernand-Mantel, A.; Mathieu, D.; Agafonov, V. Ab Initio Solid-State Heats of Formation of Molecular Salts From Ion Packing and Crystal Modelling: Application to Ammonium Crystals. Mol. Phys. 2004, 102, 253–258. (11) Perlovich, G. L.; Kurkov, S. V.; Hansen, L. K.; Bauer-Brandl, A. Thermodynamics of Sublimation, Crystal Lattice Energies, and Crystal Structures of Racemates and Enantiomers:(+)-And (±)-Ibuprofen. J. Pharm. Sci. 2004, 93, 654–666. (12) Feng, S.; Li, T. Predicting Lattice Energy of Organic Crystals by Density Functional Theory With Empirically Corrected Dispersion Energy. J. Chem. Theory Comput. 2006, 2, 149–156. (13) Maschio, L.; Civalleri, B.; Ugliengo, P.; Gavezzotti, A. Intermolecular Interaction Energies in Molecular Crystals: Comparison and Agreement of Localized Moller-Plesset 2, Dispersion-Corrected Density Functional, and Classical Empirical Two-Body Calculations. J. Phys. Chem. A 2011, 115, 11179–11186. (14) Thomas, S. P.; Spackman, P. R.; Jayatilaka, D.; Spackman, M. A. Accurate Lattice Energies for Molecular Crystals from Experimental Crystal Structures. J. Chem. Theory Comput. 2018, 14, 1614–1623. (15) Ouvrard, C.; Mitchell, J. B. Can we predict lattice energy from molecular structure? Acta Cryst. B 2003, 59, 676–685. (16) Mathieu, D. Simple Alternative to Neural Networks for Predicting Sublimation Enthalpies from Fragment Contributions. Ind. Eng. Chem. Res. 2012, 51, 2814–2819.

26

ACS Paragon Plus Environment

Page 26 of 32

Page 27 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

(17) Politzer, P.; Murray, J. S.; Grice, M. E.; Desalvo, M.; Miller, E. Calculation of heats of sublimation and solid phase heats of formation. Mol. Phys. 1997, 91, 923–928. (18) Mathieu, D.; Bougrat, P. Model equations for estimating sublimation enthalpies of organic compounds. Chem. Phys. Lett. 1999, 303, 369–375. (19) Murray, J. S.; Brinck, T.; Lane, P.; Paulsen, K.; Politzer, P. Statistically-based interaction indices derived from molecular surface electrostatic potentials: a general interaction properties function (GIPF). J. Mol. Struct. (Theochem) 1994, 307, 55–64. (20) Mathieu, D.; Simonetti, P. Evaluation of solid-state formation enthalpies for energetic materials and related compounds. Thermochim. Acta 2002, 384, 369–375. (21) Kim, C. K.; Lee, K. A.; Hyun, K. H.; Park, H. J.; Kwack, I. Y.; Kim, C. K.; Lee, H. W.; Lee, B. S. Prediction of physicochemical properties of organic molecules using van der Waals surface electrostatic potentials. J. Comput. Chem. 2004, 25, 2073–2079. (22) Byrd, E. F. C.; Rice, B. M. Improved Prediction of Heats of Formation of Energetic Materials Using Quantum Mechanical Calculations. J. Phys. Chem. A 2006, 110, 1005– 1013. (23) Anguang, H.; Brian, L.; Sergey, D.; Hakima, A.-R.; Louis-Simon, L.; Hong, G. Theoretical Prediction of Heats of Sublimation of Energetic Materials Using Pseudo-Atomic Orbital Density Functional Theory Calculations. Propellants Explos. Pyrotech. 2007, 32, 331–337. (24) Puri, S.; Chickos, J. S.; Welsh, W. J. Three-Dimensional Quantitative StructureProperty Relationship (3D-QSPR) Models for Prediction of Thermodynamic Properties of Polychlorinated Biphenyls (PCBs): Enthalpy of Sublimation. J. Chem. Inf. Comput. Sci. 2002, 42, 109–116.

27

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(25) Gharagheizi, F. A new molecular-based model for prediction of enthalpy of sublimation of pure components. Thermochim. Acta 2008, 469, 8–11. (26) Perlovich, G. L.; Raevsky, O. A. Sublimation of Molecular Crystals: Prediction of Sublimation Functions on the Basis of HYBOT Physicochemical Descriptors and Structural Clusterization. Crystal Growth & Design 2010, 10, 2707–2712. (27) Salahinejad, M.; Le, T. C.; Winkler, D. A. Capturing the Crystal: Prediction of Enthalpy of Sublimation, Crystal Lattice Energy, and Melting Points of Organic Compounds. J. Chem. Inf. Model. 2013, 53, 223–229. (28) Keshavarz, M. H.; Bashavard, B.; Goshadro, A.; Dehghan, Z.; Jafari, M. Prediction of heats of sublimation of energetic compounds using their molecular structures. J. Thermal Anal. Calorim. 2015, 120, 1941–1951. (29) Gharagheizi, F.; Sattari, M.; Tirandazi, B. Prediction of Crystal Lattice Energy Using Enthalpy of Sublimation: A Group Contribution-Based Model. Ind. Eng. Chem. Res. 2011, 50, 2482–2486. (30) Gharagheizi, F.; Ilani-Kashkouli, P.; W. E. Acree, J.; Mohammadi, A. H.; Ramjugernath, D. A group contribution model for determining the sublimation enthalpy of organic compounds at the standard reference temperature of 298 K. Fluid Phase Equil. 2013, 354, 265–285. (31) Suntsova, M. A.; Dorofeeva, O. V. Prediction of enthalpies of sublimation of highnitrogen energetic compounds: Modified Politzer model. J. Mol. Graph. Model. 2017, 72, 220 – 228. (32) Dorofeeva, O. V.; Ryzhova, O. N. Gas-Phase Enthalpies of Formation and Enthalpies of Sublimation of Amino Acids Based on Isodesmic Reaction Calculations. J. Phys. Chem. A 2014, 118, 3490–3502.

28

ACS Paragon Plus Environment

Page 28 of 32

Page 29 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

(33) Mathieu, D.; Simonetti, P. Predicting and Evaluating Performance of Energetic Salts: Models and Theoretical Tools. International J. Energetic Mater. and Chemical Propulsion 2009, 8 . (34) Gutowski, K. E.; Rogers, R. D.; Dixon, D. A. Accurate Thermochemical Properties for Energetic Materials Applications. II. Heats of Formation of Imidazolium-, 1,2,4Triazolium-, and Tetrazolium-Based Energetic Salts from Isodesmic and Lattice Energy Calculations. J. Phys. Chem. B 2007, 111, 4788–4800. (35) Glasser, L.; Jenkins, H. D. B. Volume-Based Thermodynamics: A Prescription for Its Application and Usage in Approximation and Prediction of Thermodynamic Data. J. Chem. Eng. Data 2011, 56, 874–880. (36) Glasser, L.; Jenkins, H. D. B. Predictive thermodynamics for ionic solids and liquids. Phys. Chem. Chem. Phys. 2016, 18, 21226–21240. (37) Glorian, J.; Han, K.; Braun, S.; Ciszek, F.; Baschung, B. Investigation of the Heat of Explosion of Energetic Materials: An Experimental and Computational Study. Energetic materials : synthesis, processing, performance : 49th International Annual Conference of the Fraunhofer ICT, June 26-29, 2018, Convention Center, Karlsruhe, Germany. 2018; p V15. (38) Acree, W. A.; Chickos, J. S. Phase Transition Enthalpy Measurements of Organic and Organometallic Compounds. Sublimation, Vaporization and Fusion Enthalpies From 1880 to 2015. Part 1. C1 - C10 . J. Phys. Chem. Ref. Data 2016, 45, 033101. (39) http://www.statsmodels.org. (40) Ponder, J. W. TINKER, version 7.1.2, https://dasher.wustl.edu/tinker/ (accessed March 2016).

29

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(41) Jeffrey, G. A.; Takagi, S. Hydrogen-Bond Structure in Carbohydrate Crystals. Acc. Chem. Res. 1978, 11, 264. (42) McDonagh, J. L.; Palmer, D. S.; Mourik, T. v.; Mitchell, J. B. O. Are the Sublimation Thermodynamics of Organic Molecules Predictable? Journal of Chemical Information and Modeling 2016, 56, 2162–2179. (43) Acree, W. A.; Chickos, J. S. Phase Transition Enthalpy Measurements of Organic and Organometallic Compounds. Sublimation, Vaporization and Fusion Enthalpies From 1880 to 2015. Part 2. C11 - C192 . J. Phys. Chem. Ref. Data 2017, 46, 013104. (44) Curtiss, L. A.; Redfern, P. C.; Raghavachari, K. Gaussian-4 Theory. J. Chem. Phys. 2007, 126, 084108. (45) Curtiss, L. A.; Redfern, P. C.; Raghavachari, K. Gaussian-4 Theory Using Reduced Order Perturbation Theory. J. Chem. Phys. 2007, 127, 124105. (46) Stewart, J. J. P. Optimization of parameters for semiempirical methods I. Method. J. Comput. Chem. 10, 209–220. (47) Granovsky, A. A. Firefly version 8.0. http://classic.chem.msu.su/gran/firefly/index.htm. (48) Mathieu, D.; Alaime, T. Insight into the contribution of individual functional groups to the flash point of organic compounds. J. Hazard. Mater. 2014, 267, 169–174. (49) Mathieu, D. Solubility of Organic Compounds in Octanol: Improved Predictions Based on the Geometrical Fragment Approach. Chemosphere 2017, 182, 399–405. (50) Beaucamp, S.; Mathieu, D.; Agafonov, V. Optimal partitioning of molecular properties into additive contributions: the case of crystal volumes. Acta Cryst. B 2007, 63, 277– 284. (51) Mathieu, D.; Bouteloup, R. Reliable and versatile model for the density of liquids based on additive volume increments. Ind. Eng. Chem. Res. 2016, 55, 12970–12980. 30

ACS Paragon Plus Environment

Page 30 of 32

Page 31 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

(52) Politzer, P.; Murray, J. S.; Grice, M. E.; Desalvo, M.; Miller, E. Calculation of heats of sublimation and solid phase heats of formation. Mol. Phys. 1997, 91, 923–928. (53) Rice, B. M.; Pai, S. V.; Hare, J. Predicting heats of formation of energetic materials using quantum mechanical calculations. Combus. & Flame 1999, 118, 445–458. (54) Wit, H. D.; Miltenburg, J. V.; Kruif, C. D. Thermodynamic properties of molecular organic crystals containing nitrogen, oxygen, and sulphur 1. Vapour pressures and enthalpies of sublimation. J. Chem. Thermodyn. 1983, 15, 651 – 663. (55) Sellers, P. Enthalpies of Formation of Squaric Acid and the Corresponding Diethyl Ether. Actra Chem. Scan. 1971, 25, 2194–2198. (56) Sellers, P. Enthalpy of Formation of Tetramethyl-1,3-cyclobutane-dione. Actra Chem. Scan. 1971, 25, 2291–2294. (57) Gupta, S.; Basant, N.; Singh, K. P. Three-Tier Strategy for Screening High-Energy Molecules Using StructureâĂŞProperty Relationship Modeling Approaches. Ind. Eng. Chem. Res. 2016, 55, 820–831. (58) Lai, W.-P.; Yu, T.; Liu, Y.-Z.; Ma, Y.-D.; Lian, P.; Ge, Z.-X.; Lv, J. Study on the computer-aided design of high energetic compounds based on the 1,2,3,4-tetrazine-1,3dioxide frame. J. Mol. Model. 2017, 23, 340.

31

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Graphical TOC Entry

32

ACS Paragon Plus Environment

Page 32 of 32