Accelerating the discovery of energetic melt-castable materials by a high-throughput virtual screening and experimental approach

Abstract

With the growth of chemical data, computation power and algorithms, machine learning-assisted high-throughput virtual screening (ML-assisted HTVS) is revolutionizing the research paradigm of new materials. Herein, a combined ML-assisted HTVS and experimental approach was applied to accelerate the search for energetic melt-castable materials with promising properties. The ML-assisted HTVS system is composed of high-throughput molecular generation (in a heuristic enumeration method) and five machine learning-based property prediction models (including density, melting point, decomposition temperature, detonation velocity, and detonation pressure). Using this system, we rapidly targeted 136 promising candidates from a generated molecular space containing 3892 molecules. With extensive efforts on experimental synthesis, eight new energetic melt-castable materials (MC-1 to MC-8) were obtained, and their measured properties were in good agreement with the predicted results. This work verifies the effectiveness of the combined ML-assisted HTVS and experimental approach for the accelerated discovery of energetic melt-castable materials.

---

Introduction

The demand for exploring a large number of potential new materials coupled with advances in computational power and algorithms gave birth to high-throughput virtual screening (HTVS), which features in silico high-throughput computation to help experimental chemists rapidly identify the most promising candidates for experimental endeavors.^1–3 The HTVS approach is becoming a ground-breaking tool for new materials discovery and has helped scientists accelerate the search for various new materials (such as alloy materials, metal–organic frameworks (MOFs), perovskite materials, etc.).^4–9 However, the application of HTVS in materials discovery is still in its infancy in relatively few fields, and the accelerated discovery using HTVS has not been achieved for most materials (e.g., energetic materials).

Energetic melt-castable materials are an important branch in the field of energetic materials, which exhibit a reversible solid–liquid phase transition with a relatively narrow melting range (70–120 °C)^10–12 and are typically used as the carrier to load other crystalline energetic (or non-energetic) components into the different dimensions and shapes required for applications.¹³ To date, the most widely used energetic melt-castable material is still 2,4,6-trinitrotoluene (TNT), a high-density aromatic compound first discovered 150 years ago that has many drawbacks, including a low energy, impurities, toxicity, etc.^14,15 In past decades, a few new energetic melt-castable materials, such as 1-methyl-3,4,5-trinitropyrazole (MTNP), 1,3,3-trinitroazetidine (TNAZ), and bis(1,2,4-oxadiazole)bis(methylene) dinitrate (BOM) have been reported as potential TNT replacements,^16–18 but their comprehensive performances are far from industrial applicability and are still in the laboratory evaluation stage.

In the design of new energetic melt-castable materials, four-dimensional property parameters (i.e., energy, sensitivity, melting point, and decomposition temperature) are considered simultaneously, which is more complicated than common energetic materials (which mainly concern energy and sensitivity).^19–22 In addition, there is a lack of models with good generalization to predict the melting point and decomposition temperature of energetic materials, and the molecular design of new energetic melt-castable materials primarily depends on scientific intuition, which is a significantly low-efficiency and high-uncertainty process. Therefore, it is still challenging to develop new energetic melt-castable materials with desired properties.

The development of energetic melt-castable materials can be accelerated by performing machine learning (ML)-assisted HTVS. The ML technology can provide us with statistical algorithms for training high-efficient prediction models of multi-dimensional properties^23–26 and has been widely applied in field of chemistry.^27–31 To our knowledge, ML-assisted HTVS for energetic melt-castable materials has not been established. In this work, we developed the first ML-assisted HTVS for the design of new energetic melt-castable materials, and then we synthesized the screened target molecules. Compared with traditional trial-and-error, which relies heavily on chemical intuition and iterative design-experiment feedback, the ML-assisted HTVS method is significantly more cost-effective and efficient (Fig. 1).

Fig. 1 Comparison between traditional research paradigm and ML-assisted HTVS. Typical procedures of traditional molecular discovery paradigm based on chemical intuition (left) and the accelerated discovery of new materials guided by ML-assisted HTVS (right).

Results and discussion

The framework of HTVS system

The ML-assisted HTVS system for the discovery of energetic melt-castable materials works as follows (Fig. 2). First, a molecular library with many candidates is established by high-throughput molecular generation via permutation and combination between energetic molecular skeletons and substituents (Fig. 2A and B). Then, the generated molecules are imported into ML-based property models (including density, melting point, decomposition temperature, detonation velocity, and detonation pressure) to undergo rapid and accurate property calculations. Potential molecules with high energy, low melting point, and high decomposition temperature can be filtered out by setting reasonable screening criteria (Fig. 2C). Finally, the molecules with promising comprehensive performances can be screened for rational synthesis and characterization (Fig. 2D) rather than spending a great deal of time and effort on trial-and-error.

Fig. 2 Overall framework for ML-assisted HTVS applied in this work. (A) Molecular fragments for constructing the molecular library. (B) Flow chart for generating molecular structures. (C) Schematic of ML-assisted properties prediction and molecule screening for locking target compounds with tailored properties. (D) Synthesis and characterization.

Molecular generation by heuristic enumeration method

For organic functional materials, there is a universal rule that flexible substituent groups (such as methyl and alkyl chains) can increase molecular degrees of freedom and thereby lower the melting point (Fig. S1†). Here, we propose the permutation and combination of heterocyclic molecular skeletons (five-membered azole and oxadiazole rings) and flexible substituent groups (such as methyl, methoxy groups, and nitrate/nitramine alkyl chains) to rapidly generate thousands of candidate molecules as the original molecular library. In this molecular design rule, the hydrogen atom in N–H bonds of five-membered azole rings (such as pyrazole, imidazole, triazole, and tetrazole) are all substituted by an R2 sequence of substitution groups, and the hydrogen atom in the C–H bonds of five-membered rings are partially or all substituted by an R1 and R2 sequence of substituents (Fig. 2A). This molecular generation method is named as “heuristic enumeration” because it requires initial input structures as a template and involves an enumeration process that covers all reasonable structures.^32,33

Using the pyrazole ring as an example, the flow chart of the heuristic enumeration process is shown in Fig. 2B. First, the hydrogen atom in the N–H bond of pyrazole was substituted by an R2 sequence of substituents (labelled as ⑤, ⑥, ⑦ in Fig. 2A) to generate three alkyl-substituted pyrazole derivatives. Next, one hydrogen atom (labelled by 1, 2, 3) in three C–H bonds on the alkyl-substituted pyrazole ring (Fig. 2B middle) was randomly selected and replaced by an R1 and R2 sequence of substituents (labelled as ① to ⑦ in Fig. 2A) one-by-one to enumerate all possible molecular structures (this process is defined as the first generation). After sanitization and deduplication, the first generation molecules repeated the above pick-up and substitution process to generate second-generation molecular structures with two hydrogen atoms substituted. This pick-up and substitution process can be iterated until a pre-defined substitution number is reached (usually limited to the number of hydrogen atoms in C–H bonds on five-membered rings). The whole generation process was implemented by a homemade script based on the RDkit package.³⁴ Finally, we obtained a molecular space containing 3892 total molecules (Data S1†).

Featurization, training, and validation of ML models

To train the property prediction models for the screening of promising candidates, adequate structure–property data and reasonable feature extraction are required;^35,36 however, due to high risk and special purposes, open and comprehensive data of energetic materials are still scarce. Here, we manually collected more than one thousand data from the literature to construct a basic structure–property database for ML model training. Next, molecular structures were transformed into vectors by a composite descriptor set composed by electrotopological (E-state) fingerprints and custom descriptors. The E-state fingerprints contain the number of atoms for certain types and the corresponding electrotopological state indices, which have been applied to boiling point prediction.^37,38 Custom descriptors (see Table S1†) were defined to enhance the description of molecular shapes, energetic features (such as oxygen balance and the number of nitro groups), and non-covalent interactions (such as TPSA³⁹ and partial charges). The calculations of these descriptors are mainly based on the molecular structure and forcefield method, while electronic structure calculations are deliberately avoided. The short calculation time for molecular descriptors is beneficial for the high-throughput virtual screening of molecules, even on personal computers.

During training models by kernel ridge regression (KRR) algorithm, two metrics, R² (coefficient of determination) and MAE (mean absolute error), were used to evaluate the model accuracy. The scatter plot of predictions versus observations and the deviation distribution for density, detonation velocity (D_v), melting point (T_m), and decomposition temperature (T_d) models (the detonation pressure (P) model plot was highly consistent with the D_v model, as shown in Fig. S2†) were visualized in Fig. 3A. The density and D_v models achieved good accuracy and generalization with a high R² (0.93 and 0.83) and low MAE (0.042 g cm⁻³ and 0.24 km s⁻¹). Meanwhile, in the parity plot of density and D_v, data points are evenly distributed on both sides around the diagonal, and the error distribution of the test set is consistent with that of the training set; however, the metrics for T_m (R²: 0.60 and MAE: 32.5 °C) and T_d (R²: 0.62 and MAE: 30.8 °C) are less satisfactory. Even so, melting point prediction is still complicated, and our T_m model reached a high accuracy in comparation with former reports (Fig. S3†). In addition, T_d models for organic molecules are seldom investigated.⁴⁰ Overall, our trained T_m and T_d models provide a powerful tool for predicting the thermal properties of energetic materials. The reason for the relatively low accuracy of T_m and T_d models is probably assigned to that our descriptors are still relatively insufficient for describing the intermolecular interactions and the complicated physical picture during melting or decomposition, as well as the experimental & chemical uncertainties.⁴¹ Fortunately, the MAE for the T_m and T_d models are very close; hence, we used the difference between T_m and T_d (temperature range, T_r) as an additional filter condition. We wish the effect of error cancellation will make screening results more instructive.

Fig. 3 Evaluation on the trained models and feature importance for T_m model. (A) Parity plot of prediction values versus observation values and error distribution for different models (all trained by kernel ridge regression algorithm, KRR). Purple and green colors indicate the training and test set, respectively. (B) Evaluation on feature importance for T_m models from random forest (RF) algorithm.

To determine which descriptor has vital influence on the melting point, we applied a random forest (RF) algorithm to train the T_m model and return the feature importance for custom descriptors (Fig. 3B).⁴² The RF model metrics are 0.61 (R²) and 31.0 °C (MAE), respectively, which are close to those of the KRR model. From the feature importance, we can clearly see that descriptors related to intermolecular interactions, such as TPSA (topological polar surface area, 10.34%), nR (number of rings, 6.64%), Min/Max partial charge (3.81%/3.41%), PBF (the plane of best fit, 3.53%), and nH bond (number of hydrogen bonds, 3.49%), are more important for predicting T_m. This result corroborates the argument that intermolecular interactions have an obvious effect on the melting point. In addition, the sum of the importance of custom descriptors (58.38%) exceeds the proportion of E-state fingerprints, which indicates our custom descriptors are also essential for training the T_m model.

Property prediction and screening of generated molecules

With well-trained property prediction models in hand, we attempted to predict the properties of generated molecules (a total of 3892 five-membered nitrogen heterocyclic molecules) and transfer the chemical space to the property space. The initial property space is visualized by a pair–grid plot (Fig. S4†) and violin plots (Fig. 4A) which are categorized by heterocyclic molecular skeletons. In the pair–grid plot, we can see that density shows a strong positive correlation with D_v (R = 0.97, R is short for correlation coefficient) and P (R = 0.99), while T_d shows a negative correlation with density (R = −0.63), D_v (R = −0.59), and P (R = −0.64). The results are in accordance with the commonly observed phenomenon that density can affect the energy of energetic materials, but thermostability usually decreases as the energy increases,^43,44 which implies the effectiveness of our ML models. Meanwhile, T_m showed weakly positive correlations with energy properties including density (R = 0.46), D_v (R = 0.42), and P (R = 0.49). This is understandable because a high energy (especially at a high density) indicates the close molecular packing in crystals with a high possibility of forming stronger intermolecular interactions, which always increase the T_m. Contrary to our expectations, there was no obvious correlation between T_m and T_d (R = −0.04). This phenomenon may indicate that T_m and T_d are totally determined different physical factors (or features from the view of ML, Fig. S5†), which makes it possible for us to identify compounds with desirable melting point range (70–120 °C) and large T_r.

Fig. 4 Screening and synthesis of molecules. (A) Violin plot for all generated molecules categorized by heterocyclic molecular skeletons. White points, heavy lines, thin lines, and colored areas indicate the median, interquartile range, 95% confidence interval, and probability density, respectively. (B) Screening process; the table records the number of molecules in each screening step. (C) Eight synthesized molecules. (D) Schematic of synthetic routes.

The violin plots show that molecules with different heterocyclic molecular skeletons exhibit various property preferences (Fig. 4A). For azole ring-based molecules, their energy properties (such as density, D_v and P) and thermostability (T_d) show positive and negative relationships with the numbers of nitrogen atoms in the azole rings, respectively. For example, pyrazole and imidazole-based molecules tend to be thermally stable, while their energies (D_v and P) are also relatively low. In contrast, tetrazole-based molecules usually exhibit high energies (D_v and P) but the lowest T_d. These results support the universal rule that a high nitrogen content in the molecular skeleton increases the energy but decreases the thermal stability.⁴⁵ Due to the introduction of oxygen atoms in energetic parent rings, oxadiazole-based molecules usually have a high density and energy (D_v and P); however, their melting points are relatively low, with the median around 60 °C (Oxa-1 to Oxa-4 in Fig. 4A), which are lower than the desired melting point (70–120 °C) for energetic melt-cast materials. From the view of melting points, pyrazole and imidazole-based molecules are expected to be the most promising candidates because their median of predicted melting points are around 85 °C (Pyr and Imi in Fig. 4A).

To rapidly and precisely identify the most promising energetic melt-castable materials, we narrowed the property space by introducing additional criteria (density >1.65 g cm⁻³, T_m <110 °C, T_d >190 °C, and T_r >90 °C) step-by-step. The number and degree of change for each screening step are summarized in Fig. 4B. The initial molecular library contains 3892 molecules, and molecules derived from pyrazole and imidazole accounted for around 79% (1533 and 1533) of the total amount because both of these parent rings possess four substitutable sites for the construction of new molecules. After a four-step screening process, the number of satisfied molecules shrank to 136 (almost 29 times smaller than the initial molecular library), and the remaining molecules derived from pyrazole and imidazole still occupied the majority (total 89.7%). Although the D_v and P criteria were not explicitly defined during the screening process, the D_v and P of 136 screened molecules separately exceeded 7500 km s⁻¹ and 23 GPa, which are superior to those of TNT (D_v: 7.067 km s⁻¹ and P: 19.3 GPa calculated by our ML models). These results indicate that the set screening criteria are reasonable and can efficiently target potential molecules with good comprehensive performance properties (including a high energy, low melting point, and high decomposition temperature). Details about the properties of screened molecules are summarized in Table S2.†

Experimental validation and property studies

Although 136 potential molecules were screened, we think around twelve molecules are valuable for synthetic endeavor after eliminating structural isomers, reported molecules, and synthesis-inaccessible ones (Fig. S6†). Significant efforts were made to synthesize these molecules, and eight new energetic melt-castable molecules (MC-1 to MC-8 in Fig. 4C) were obtained through reasonable reaction routes (Fig. 4D). Detailed synthesis descriptions and structural characterizations can be seen in ESI.†

As expected, the eight compounds all exhibited the desired melting endothermic peaks below 110 °C, although the ones of MC-2 (T_m: 68.4 °C) and MC-3 (T_m: 65.2 °C) are slightly lower than 70 °C (Fig. 5). For the development of energetic melt-castable materials, a high decomposition temperature and large T_r values are the most important thermal properties besides a reasonable melting point range (70–120 °C). Among the eight energetic compounds, the thermostability of MC-1, MC-3, and MC-5 are good, with onset decomposition temperatures exceeding 170 °C and T_r values all higher than 100 °C. Combined with the fact that the melting point (65.2 °C) of MC-3 is lower than 70 °C (the lower limit of the reasonable melting point range), we think that MC-1 and MC-5 have great potential as energetic melt-castable materials due to their thermal properties. The MAE of the melting point and decomposition temperature in the ML models were 32.5 °C and 30.8 °C, respectively (Fig. 3A), but the actual deviations between the experimental and predicted values of these eight compounds are all below 30 °C, and most are lower than 20 °C (Fig. 6b and c). These results clearly indicate that our ML models can generally predict the T_m and T_d of novel compounds using our training data.

Fig. 5 Thermal behavior and crystal structures of synthesized compounds. DSC/TG curves for synthesized compounds at heating rate of 10 °C min⁻¹; insets are the respective molecular (left inset) and crystal structures (right inset). The tangent intersection is the decomposition onset.

Fig. 6 Comparison between predicted and standard values. (A) Densities, (B) melting points, (C) decomposition peak temperatures, (D) detonation velocities, and (E) detonation pressures of eight melt-castable materials. (F) Comprehensive property comparison of MC-1, MC-5 and TNT, all values are normalized to those of TNT.

After slow solvent evaporation at room temperature, single crystals of these eight energetic compounds were obtained (inset pictures in Fig. 5, S7 and Tables S3–S10†), and their crystal densities were in the range of 1.68–1.83 g cm⁻³ (red dots in Fig. 6A), which are higher than that of TNT (1.65 g cm⁻³). Their detonation velocities and detonation pressures (evaluated by Explo5 (v6.02) program) ranged from 7.843 km s⁻¹ to 8.705 km s⁻¹ and 25.2 GPa to 33.8 GPa (red dots in Fig. 6D and E),⁴⁶ respectively, which are also both higher than those of TNT (7.303 km s⁻¹ and 21.3 GPa). Overall, these eight compounds exhibited much higher energies than TNT, suggesting they can be used as high-energy melt-castable materials. More importantly, the densities and detonation performances predicted by our ML models are very similar to the experimental and calculated values (Fig. 6a, d and e), demonstrating the ability of our ML models to predict the energy properties of new energetic materials. Except for MC-1 (IS = 8 J), which contains a nitramine group in its molecular structure, the remaining seven compounds exhibited lower impact sensitivities (IS > 20 J) than TNT (15 J). However, the friction sensitivities (ranging from 96 N to 160 N) of the eight new compounds are much higher than that of TNT (353 N) (Table S11,† details about the physicochemical properties are also summarized.), which may be attributed to their mix-crossed crystal packing structures and the lack of strong hydrogen-bond interactions (Fig. S7†).⁴⁷ After a detailed evaluation of the comprehensive properties of these eight compounds, we think that MC-1 and MC-5 show promising potential as TNT alternatives due to their reasonable melting points (75.9 °C and 88.2 °C, respectively), high decomposition temperatures (208.7 °C and 219.8 °C, respectively), large differences between melt point and decomposition peak temperature (>100 °C), and desirable energy (8.139 km s⁻¹, 28.1 GPa and 7.912 km s⁻¹, 26.3 GPa, respectively) properties (Fig. 6F).

Conclusions

A ML-assisted HTVS system was developed and applied to accelerate the search for new energetic melt-castable materials. In this HTVS system, a heuristic enumeration method and five well-trained ML models (including density, detonation velocity, detonation pressure, melting point, and decomposition temperature) were responsible for the high-throughput molecular generation and fast property prediction of energetic melt-castable materials, respectively. Using this HTVS system, we identified eight new energetic melt-castable compounds from a large chemical space (containing 3892 molecules) with high efficiency and their measured properties are in good agreement with the predicted results. Especially for compounds MC-1 and MC-5, they exhibit good comprehensive performances (including high energy, reasonable melting point and good thermostability) and are thought as good potential candidates for replacing TNT. This work demonstrates a promising research paradigm for the accelerated discovery of new energetic melt-castable materials with desired properties.

Methods

The calculation of features (molecular descriptors) including custom descriptors and electrotopological fingerprints were accomplished by RDkit library. Property models were trained by kernel ridge regression algorithm implemented in Scikit-learn package.⁴⁸ The generation of molecules was achieved by home-made scripts based on the function of RDkit. More details can be found in ESI.† The preparation and characterization of compounds can also be found in ESI.†

Author contributions

Q. Z. and Y. W. designed and supervised the study. S. S. and F. C. collected data and write the related codes. S. S., Y. W., and F. C. performed the reactions, measurements and data analysis. M. Y. and K. C. W. performed the single-crystal measurement and crystallographic structural analysis. S. S., Y. W. and Q. Z. prepared the manuscript.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We thank Science Challenge Project (TZ2018004) and National Natural Science Foundation of China (No. 22075259, 21875228) for financial support.

Notes and references

1. E. O. Pyzer-Knapp, C. Suh, R. Gómez-Bombarelli, J. Aguilera-Iparraguirre and A. Aspuru-Guzik, Annu. Rev. Mater. Res., 2015, 45, 195–216.
2. S. Hirata and K. Shizu, Nat. Mater., 2016, 15, 1056–1057.
3. J. Noh, G. Ho Gu, S. Kim and Y. Jung, Chem. Sci., 2020, 11, 4871–4881.
4. F. Ren, L. Ward, T. Williams, K. J. Laws, C. Wolverton, J. Hattrick-Simpers and A. Mehta, Sci. Adv., 2018, 4, eaaq1566.
5. M. Fernandez, P. G. Boyd, T. D. Daff, M. Z. Aghaji and T. K. Woo, J. Phys. Chem. Lett., 2014, 5, 3056–3060.
6. Z. Li, Q. Xu, Q. Sun, Z. Hou and W. J. Yin, Adv. Funct. Mater., 2019, 29, 1807280.
7. R. Gómez-Bombarelli, J. Aguilera-Iparraguirre, T. D. Hirzel, D. Duvenaud, D. Maclaurin, M. A. Blood-Forsythe, H. S. Chae, M. Einzinger, D. G. Ha, T. Wu, G. Markopoulos, S. Jeon, H. Kang, H. Miyazaki, M. Numata, S. Kim, W. Huang, S. I. Hong, M. Baldo, R. P. Adams and A. Aspuru-Guzik, Nat. Mater., 2016, 15, 1120–1127.
8. H. Ruan, C. Yu, X. Niu, W. Zhang, H. Liu, L. Chen, R. Xiong, Q. Sun, C. Jin, Y. Liu and L. Lai, Chem. Sci., 2021, 12, 3004–3016.
9. A. Vriza, A. B. Canaj, R. Vismara, L. J. K. Cook, T. D. Manning, M. W. Gaultois, P. A. Wood, V. Kurlin, N. Berry, M. S. Dyer and M. J. Rosseinsky, Chem. Sci., 2021, 12, 1702–1719.
10. Q. Ma, Z. Zhang, W. Yang, W. Li, J. Ju and G. Fan, Energetic Materials Frontiers, 2021, 2, 69–85.
11. L. A. Wingard, P. E. Guzmán, E. C. Johnson, J. J. Sabatini, G. W. Drake and E. F. Byrd, ChemPlusChem, 2017, 82, 195–198.
12. J. C. Bennion, P. G. Lafond and J. A. Ciezak-Jenkins, Propellants, Explos., Pyrotech., 2019, 44, 1015–1020.
13. F. Chen, S. Song, Y. Wang, Y. Liu and Q. Zhang, Energetic Materials Frontiers, 2020, 1, 157–164.
14. M. Anniyappan, V. K. Vijay, R. S. Amit and J. K. Nair, J. Energ. Mater., 2020, 38, 111–125.
15. F. Abrishami, M. Chizari, N. Zohari and S. A. Pourmosavi, Propellants, Explos., Pyrotech., 2019, 44, 1446–1449.
16. G. Hervé, C. Roussel and H. Graindorge, Angew. Chem., Int. Ed., 2010, 49, 3177–3181.
17. N. Sikder, A. K. Sikder, N. R. Bulakh and B. Gandhe, J. Hazard. Mater., 2004, 113, 35–43.
18. E. C. Johnson, J. J. Sabatini, D. E. Chavez, R. C. Sausa, E. F. Byrd, L. A. Wingard and P. E. Guzmàn, Org. Process Res. Dev., 2018, 22, 736–740.
19. J. Zhang, L. A. Mitchell, D. A. Parrish and J. M. Shreeve, J. Am. Chem. Soc., 2015, 137, 10532–10535.
20. D. E. Chavez, M. A. Hiskey and R. D. Gilardi, Angew. Chem., Int. Ed., 2000, 39, 1791–1793.
21. M. C. Schulze, B. L. Scott and D. E. Chavez, J. Mater. Chem. A, 2015, 3, 17963–17965.
22. V. W. Manner, M. J. Cawkwell, E. M. Kober, T. W. Myers, G. W. Brown, H. Tian, C. J. Snyder, R. Perriot and D. N. Prestona, Chem. Sci., 2018, 9, 3649–3663.
23. K. T. Butler, D. W. Davies, H. Cartwright, O. Isayev and A. Walsh, Nature, 2018, 559, 547–555.
24. J. B. O. Mitchell, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2014, 4, 468–481.
25. D. C. Elton, Z. Boukouvalas, M. S. Butrico, M. D. Fuge and P. W. Chung, Sci. Rep., 2018, 8, 9059.
26. P. Kang, Z. Liu, H. Abou-Rachid and H. Guo, J. Phys. Chem. A, 2020, 124, 5341–5351.
27. M. H. S. Segler, M. Preuss and M. P. Waller, Nature, 2018, 555, 604–610.
28. A. D. Casey, S. F. Son, I. Bilionis and B. C. Barnes, J. Chem. Inf. Model., 2020, 60, 4457–4473.
29. F. Brockherde, L. Vogt, L. Li, M. E. Tuckerman, K. Burke and K. R. Müller, Nat. Commun., 2017, 8, 872.
30. K. Fujimura, A. Seko, Y. Koyama, A. Kuwabara, I. Kishida, K. Shitara, C. A. J. Fisher, H. Moriwake and I. Tanaka, Adv. Energy Mater., 2013, 3, 980–985.
31. M. K. Nielsen, D. T. Ahneman, O. Riera and A. G. Doyle, J. Am. Chem. Soc., 2018, 140, 5004–5008.
32. R. Gani and E. A. Brignole, Fluid Phase Equilib., 1983, 13, 331–340.
33. M. Sumita, X. Yang, S. Ishihara, R. Tamura and K. Tsuda, ACS Cent. Sci., 2018, 4, 1126–1133.
34. G. Landrum, RDKit: Open-source cheminformatics, 2016.
35. P. Pankajakshan, S. Sanyal, O. E. de Noord, I. Bhattacharya, A. Bhattacharyya and U. Waghmare, Chem. Mater., 2017, 29, 4190–4201.
36. J. Townsend, Nat. Commun., 2020, 11, 3230.
37. L. H. Hall and L. B. Kier, J. Chem. Inf. Comput. Sci., 1995, 35, 1039–1045.
38. L. H. Hall and C. T. Story, J. Chem. Inf. Comput. Sci., 1996, 36, 1004–1014.
39. S. Prasanna and R. J. Doerksen, Curr. Med. Chem., 2009, 16, 21–41.
40. Q. Wang, J. Wang and M. D. Larranaga, J. Therm. Anal. Calorim., 2013, 111, 1033–1037.
41. G. Sivaraman, N. E. Jackson, B. Sanchez-Lengeling, Á. Vázquez-Mayagoitia, A. Aspuru-Guzik, V. Vishwanath and J. J. de Pablo, Chemrxiv, 2019, 9914378.v1.
42. W. Hu, S. Ye, Y. Zhang, T. Li, G. Zhang, Y. Luo, S. Mukamel and J. Jiang, J. Phys. Chem. Lett., 2019, 10, 6026–6031.
43. X. X. Zhao, S. H. Li, Y. Wang, Y. C. Li, F. Q. Zhao and S. P. Pang, J. Mater. Chem. A, 2016, 4, 5495–5504.
44. H. Li, L. Zhang, N. Petrutik, K. Wang, Q. Ma, D. Shem-Tov, F. Zhao and M. Gozin, ACS Cent. Sci., 2019, 6, 54–75.
45. V. P. Sinditskii, V. Y. Egorshev, G. F. Rudakov, A. V. Burzhava, S. A. Filatov and L. D. Sang, Thermochim. Acta, 2012, 535, 48–57.
46. M. Sućeska, EXPLO5–Computer program for calculation of detonation parameters, in Proc. of 32nd Int. Annual Conference of ICT, Karlsruhe, Germany, 2001.
47. C. Zhang, F. Jiao and H. Li, Cryst. Growth Des., 2018, 18, 5713–5726.
48. F. Pedregosa, A. G. Gaël Varoquaux, M. Vincent, T. Bertrand, G. Olivier, B. Mathieu, P. Peter, W. Ron and D. Vincent, J. Mach. Learn. Res., 2011, 12, 2825–2830.

---

Footnotes

† Electronic supplementary information (ESI) available. CCDC 2047538, 2047548, 2047557, 2047549, 2047558, 2047559, 2056920 and 2047550. For ESI and crystallographic data in CIF or other electronic format see DOI: 10.1039/d1ta04441a

‡ These authors contributed equally to this work.

---