Sequential construction of stable nitrogen–oxygen compounds using high-throughput quantum mechanical calculations and customized machine learning model

Abstract

Nitrogen and oxygen are the two most abundant elements in the atmosphere, yet stable compounds composed solely of these elements are relatively scarce. Conceiving novel stable nitrogen–oxygen compounds remains a formidable challenge for current experimental and theoretical research. In this study, we developed a sequential construction strategy to design 168 nitrogen–oxygen compounds with distinct structural innovation, followed by high-throughput quantum mechanical calculations with the highest possible accuracy. From the resulting 7820 structural and property parameters, we created a customized machine learning model that outperforms universal models in accuracy with 13.8% greater robustness across various data splits, achieving stable and high performance on small datasets. Data-driven analysis revealed the energy and electron-related characteristics as key factors in regulating thermodynamic stability, while physics-driven insights uncovered that electron delocalization and hyperstatic constraints fine-tune mechanical firmness. Among the designed nitrogen–oxygen compounds, 106 are expected to be more stable than the known compound N₂O₄, out of which 61 are expected to be even more stable than N₂O₅. Furthermore, their energy densities surpass those of all currently used nitrogen–oxygen oxidizers by 8.3–16.8%, highlighting our newly proposed compounds potential for use in rocket bipropellant systems. Our developed machine learning platform features a user-friendly graphical interface for easy assessment and may be of interest to researchers in other fields, including chemical industry and energy sectors.

---

1. Introduction

Nitrogen and oxygen are the two most abundant elements in our atmosphere, comprising approximately 75.5% and 23.1% of its total mass, respectively.¹ Despite their prevalence, stable compounds comprised solely of these two elements are relatively scarce. To date, only 12 pure oxygen and neutral nitrogen–oxygen compounds have been detected or synthesized, accounting for less than 0.01% of other inorganic oxidizing agents.^2–13 Notable examples include compounds of NO,^14,15 N₂O,^16,17 N₂O₄,^18,19 and N₂O₅,^20,21 wherein N₂O is known as laughing gas. Other species mainly exist in aqueous solutions of nitric acid (HNO₃),²² nitrous acid (HNO₂),¹⁹ (H₂N₂O₂),²³etc. Novel stable neutral nitrogen–oxygen compounds have the potential to facilitate the eco-friendly manufacturing of explosives for construction and mining, and act as high-performing oxidizers in propellants for super-long-flight space exploration. Once developed, these compounds could drive advancements in industry, construction, and aerospace science, significantly contributing to human progress.

In the context of energetic oxidizers, compounds such as N₂O and N₂O₄ are widely used as oxidizers in rocket bipropellant systems, as they not only withstand relatively high temperatures before decomposition, but also release substantial energy when extracting electrons from fuels. On the one hand, dense oxygen enhances the electron-extracting capability, increasing energy release during exothermic reactions and generating greater external work, which can result in longer flight ranges at a given mass. For example, in principle, O₄ doubles the specific impulse (I_sp) in comparison with liquid oxygen (LOX, O₂), reaching 500–700 seconds.²⁴ On the other hand, nitrogen addition capitalizes on the significant energy difference between single/double N–N bonds in the reactants and the triple N≡N bonds in the products, further enhancing energy output.²⁵ Recently, the thermostability of dinitramide (DN, –N(NO₂)₂) has been significantly improved, with a decomposition temperature (T_d) of 243 °C.²⁶ The successful synthesis of nitrogen-rich oxidizer trinitramide (TNA, N(NO₂)₃) can theoretically allow achieving the highest density impulse ever, 33% higher than that of the LOX.²⁷ Therefore, conceiving novel stable nitrogen–oxygen compounds is a highly promising and innovative research direction oriented towards the introduction of conceptually new oxidizers for high-performance aerospace applications.

However, conceiving such novel stable nitrogen–oxygen compounds poses a formidable challenge for the current state-of-the-art experimental and theoretical research. Experimentally, the intrinsic nature of nitrogen and oxygen favors their separation, while the retro-synthetic analysis of nitrogen–oxygen compounds is extremely elusive. In comparison to the separated forms of dinitrogen (N₂) and dioxygen (O₂), nitrogen–oxygen compounds may exhibit significantly higher energy levels, due to induced strain, with N–N single bond energy being 50% higher, and O–O bond energy 41% higher, posing the thermodynamic difficulties. Achieving synthesis of such nitrogen–oxygen compounds may require special reaction conditions, such as high compression and laser/microwave irradiation, to provide the activation energy necessary to overcome barriers of approximately 50–70 kJ·mol⁻¹,²⁸ indicating kinetic challenges.²⁴ It is possible that straightforward and systematic synthetic methodologies may not lead to the convenient preparation of the proposed nitrogen–oxygen compounds and “out-of-the-box” thinking is required. Theoretically, current computational algorithms for virtual compound generation exhibit evident inheritance characteristics, but the scarcity of available nitrogen–oxygen building blocks severely limits the diversity of combinations. Among the few virtually generated structures, unstable compounds vastly outnumber stable ones, further complicating the challenge to discover novel stable nitrogen–oxygen compounds.

Recent developments in high-throughput quantum mechanical calculations and machine learning methods, which have greatly enhanced the development of pharmaceuticals, catalysts, batteries, and other domains,^29–34 show clear advantages in addressing our problem. These methods' application in structural design, performance prediction, and prospective candidate screening has accelerated the synthesis of new advanced energetic materials.^35–42 However, universal machine learning models present limitations, when applied to these specialized materials, exhibiting fluctuating accuracy across different data splits.^43–45 This highlights the need for advanced machine learning models specifically tailored to predict the stability of energetic materials.

In this work, we developed a sequential construction strategy to design 168 nitrogen–oxygen compounds, free from the constraint of structural analogy, which significantly enhances structural innovation. Subsequently, we conducted high-throughput quantum mechanical calculations with the highest possible accuracy. Utilizing the resulting 7820 structural and property parameters, we created a customized machine-learning model that demonstrated superior accuracy and robustness in predicting the stability of nitrogen–oxygen compounds. Our findings reveal that energy and electron-related characteristics are key in regulating thermodynamic stability, while electron delocalization and hyperstatic constraints fine-tune mechanical firmness. Among the designed nitrogen–oxygen compounds, 106 are expected to be more stable than N₂O₄, and 61 could be even more stable than N₂O₅. Furthermore, the energy densities surpass those of all currently used oxidizers by 8.3–16.8%, highlighting their potential for use in rocket propulsion.

2. Methods

2.1. Sequential construction of nitrogen–oxygen compound structures

We propose a sequential construction strategy to design as comprehensive coverage of possible nitrogen–oxygen compound structures. Our first objective was the design of unit block-planar and cyclic fragments, ranging three-, four-, five-, and six-membered rings. Specifically, monocycles were systematically enumerated using combinatorial approaches, with nitrogen and oxygen atoms positioned in all possible positions within 3-, 4-, 5-, and 6-membered rings. Certain cyclic structures possess additional strain energy, which has the potential to enhance energy density, while cyclic planar conformations can enlarge electron delocalization, promoting thermodynamic stability. The second objective of this study was the sequential arrangement of the unit blocks, with direct and bridge connections, such as –N=N–, –O–, and –O–O– hinges, delivering backbones with various dimensions and symmetries. Higher dimensions may generally lead to increased strain energy and improved energy density, while higher symmetry helps balance atomic stress, enhancing structural stability. The third design element was the framework decoration, with one or more functional groups traversing in –N₃ and –N(NO₂)₂ functional groups. These groups feature large and positive enthalpy of formation, which are 87 and 185 kJ·mol⁻¹, respectively, further elevating the energy density of the proposed compounds. Due to virtually infinite possibilities for combination and decoration, exhaustive enumeration was not feasible. However, a search of SciFinder database revealed that nitrogen and oxygen-composed molecules with more than 35 atoms were rarely synthesized. Therefore, we focused our research on molecules containing fewer number of atoms.

Using the proposed sequential construction strategy, 168 nitrogen–oxygen compound structures were constructed (Fig. 1). With increased complexity, the three-, four-, five-, and six-membered cyclic fragments were progressively incorporated as building blocks. Each type of cyclic basic building blocks is connected either directly or by a bridge to other functional moieties, sequentially forming singular, fused, or linked skeletons, which exhibit more advanced arrangements. It is noted that acyclic compounds are constructed by direct combinations of functional groups with non-cyclic building blocks. The construction strictly adheres to chemical bonding principles, with three or five bonds for nitrogen and two for oxygen.

Fig. 1 Sequential construction of 168 nitrogen–oxygen compound structures with increasing complexity. (a) Acyclic compounds, (b) three-, (c) four-, (d) five-, and (e) six-membered cyclic primitives and their sequential constructions. 108 compounds in yellow-filled frames represent calculated stable ones, and the residual 48 in gray-filled frames are unstable. 12 compounds in purple-filled frames illustrate calculated stable structures, which have been experimentally detected or synthesized.

We note that although virtual molecule generation is particularly popular in literature, such computational algorithm-dependent approaches exhibit obvious inheritance characteristics. The innovative molecular structures usually account for the minority among the generated molecules, while the unstable molecule structures may take a vast majority. In contrast, the sequential construction strategy we proposed is free from the constraints of structural analogy and undoubtedly demonstrates significant structural innovation. One of our guiding principles is to assemble molecules from commonly stable cyclic building blocks with as few atoms as possible, connected either by thermodynamically favorable bridged linkers or through atomic fusion. Furthermore, the manual construction of the molecules is based on a comprehensive assessment of potentially stable structures, as shown by the quantum mechanical calculations in the later section. This design philosophy aimed to prioritize chemically sound structures with higher symmetry and potential stability, while minimizing the number of compounds required for high-accuracy quantum mechanical calculations. Moving forward, targeted generative models capable of exploring the full design space will significantly advance molecular design.

2.2. High-throughput quantum mechanical calculations

All the conceived compound structures were optimized using quantum mechanics methods at the B3LYP accuracy level with 6–311++G(d,p) basis set using the Gaussian 09 package.⁴⁶ Structures that fragmented upon optimization were classified as unstable and marked in gray in Fig. 1. Conversely, those that remained intact were considered potentially stable and are marked in yellow. The energetics calculations were conducted with CCSD(T) with jun-cc-pvtz basis set for most of the molecules. For other molecules that are too large to handle, the energy calculations were carried out using the MP2/6-311++G(d,p) level. A comparative calculation employing the lower-level B3LYP method revealed a substantial discrepancy ranging from 3.1 to 25.7 eV across different compounds (Fig. S4; ESI†). Notably, such significant energy deviations can critically impact the evaluation of thermodynamic stability. To ensure reliable predictions, energy calculations should be carried out with the highest attainable accuracy, and low-accuracy level methods should be avoided.

To evaluate the factors affecting the stability of the compounds, a high-throughput calculation flow was built to obtain energetic, electronic, bonding, and structural properties. The energy-related property was characterized by their energy gap, namely, the difference between the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO). The electronic property was measured by the dipole moment and the ratio of the positive average value to the negative average value (V_s⁺/V_s⁻). The bonding properties include Laplacian bond order for each chemical bond and the proportion of the O-connected bond. Laplacian bond order was calculated by integrating the Laplace function of the electron density in the bonding region using the Multiwfn package.⁴⁷ The structure-related property was evaluated via bond lengths and bond angles between the constituent atoms for each of the studied compounds, bond types and numbers, as well as the corresponding compound size, namely, the distances from the molecule's centroid to the outermost oxygen atoms.

2.3. Customized machine learning model

Machine learning model training and testing were conducted on the Alink platform developed by Alibaba Group.⁴⁸ Following the abovementioned methods, the high-throughput calculation generated 7820 property parameters of the 168 nitrogen–oxygen compounds-compassing not only energy gap, dipole moment, Laplacian bond order, and O-connected bond proportion, but also surface electrostatic potential, weakest bond length, and so on. The main parameters are provided in the ESI Spreadsheet.† These data were classified into parameters affiliated with the 120 stable compounds and those 48 unstable ones. The primary factors in regulating the target property of the nitrogen–oxygen compounds were extracted via feature importance analysis.

We developed a customized machine learning model to address the known accuracy problem associated with predicting the stability of energetic materials, particularly with small datasets. The main concept of the customized model is to evaluate the most popular universal models in dealing with classifying issues and weighted mix the best-performing ones, improving adaptability to nitrogen–oxygen compound datasets. The considered universal models include factorization machines (FM), logistic regression (LR), naive bayes (NB), support vector machine (SVM), random forest (RF), and gradient boosting decision trees (GBDT). Two main factors guided the model's development: prediction accuracy and robustness across different data splitting.

From the accuracy perspective, grid search was employed to optimize hyperparameters to improve the metrics, including Accuracy, Recall, and F1 score, as defined by

and

where in TP, FP, TN, and FN are short for true positives, false positives, true negatives, and false negatives, respectively, and Precision is the proportion of TP among all positives.

From the robustness perspective, the accuracy of the customized model was evaluated against different data splitting and compared with the aforementioned six universal models. The training and testing data were divided 100 times in different ways, and the accuracy metrics were recorded at each cycle of model training, wherein the 5-fold cross-validation method was used to mitigate the overfitting or underfitting risk. Among the customized and six universal models, the one exhibiting the minimum difference between the corresponding highest and lowest metrics demonstrates the best robustness against data splitting, implying the highest adaptation to the small dataset of nitrogen–oxygen compounds. The hyper-parameters of the customized model are detailed in Table S1 in ESI.†

2.4. Physics-driven feature selection and target-property determination

Feature selection is critical for accurate prediction of nitrogen–oxygen compounds stability, not only delivering higher accuracy metrics, but also delineating reliable structure-stability relationships. The chosen features encompass energetic, electronic, and bonding characteristics, and the values were obtained from our high-throughput calculations of the nitrogen–oxygen compounds.

Preliminary exploration was conducted to choose different combinations of feature descriptors. The results suggested that the four feature descriptors of energy gap, dipole moment, Laplacian bond order, O-connected bond proportion present the best availability and relevance to the target property, generating the highest accuracy metrics in predicting nitrogen–oxygen compound's stability.

The target property of a designed structure is its stability, assessed from both thermodynamic and mechanical perspectives. From a thermodynamic perspective, stability is assessed by structural optimization, which is currently one of the most prevalent and reliable approaches.^49–53 Compounds in which chemical bonds were fragmented upon structural optimization were classified as unstable, while those that remained intact were considered potentially stable. Compounds predicted to be thermodynamically stable could be considered feasible for synthesis. From a mechanical perspective, stability has been demonstrated to closely related to the firmness of the structure's individual bond,^54,55 which can be quantitatively characterized by the Laplacian bond order (LBO).^56,57 A compound is considered to possess mechanical stability if its minimum Laplacian bond order exceeds the threshold of 0.07. Mechanically stable compounds are considered as likely to be synthetically accessible, enabling potential applications.⁵⁸

With both thermodynamic and mechanical stability meeting application requirements, the higher stored energy density enables greater external work generation and may potentially serve as energetic oxidizers in propellants for super-long-flight space exploration. The energy density was characterized by specific impulse and heat of formation. The specific impulse was calculated using EXPLO5 v6.05 package, defined as ⁵⁹

I_sp ≈ [kT_c(1−P_aM/P_c)/M]^0.5, (4)

where T_c is the combustion chamber temperature, P_a and P_c are pressures in the combustion chamber and at the exhaust nozzle exit, respectively, M is the average molecular mass of the gaseous combustion products, and the coefficient k characterizes their heat-transfer properties.⁶⁰ The heat of formation was calculated using an atomization method at the G2 accuracy level with a sublimation correction, which was proposed by Politzer and Rice.⁶¹

3. Results

3.1 Accuracy verification of calculation: comparison to experimental data

According to the quantum mechanics calculation results, 70% of the conceived structures were thermodynamically stable, demonstrating a high probability of the potential existence of these novel nitrogen–oxygen compounds. Bond lengths and angles were measured and listed in the ESI Spreadsheet† for each compound.

Notably, 12 of our designed compounds have been previously known or experimentally synthesized,^2–13 namely, O₂, NO, O₃, NO₂, N₂O, N₂O₂, N₂O₃, N₂O₄, N₂O₃, N₄O₂, N₄O, and N₂O₅. Their structures have been characterized by X-ray diffraction and 55 bonding values were reported, as illustrated by purple-filled frames in Fig. 1. When compared to the experimental results, our calculations yield a bond length correlation coefficient (R²) of 0.98 and a root mean square error (RMSE) of 0.04 Å. The bond angle R² is 0.99 with an RMSE of 1.4°, as shown in Fig. 2(a). Additionally, 11 of our designed compounds have been previously calculated in other studies, with 56 corresponding bonding parameters reported in the literature.^2–11 However, these prior studies employed somehow lower accuracy levels. Our comparative analysis shows an R² value of 0.99 and an RMSE of 0.02 Å for bond length statistics, as well as an R² value of 0.99 and an RMSE of 0.75° for bond angle statistics, as presented in Fig. 2(b). The minimal discrepancies between our calculation outcomes and the previously reported experimental/computational values strongly support our computational method's high accuracy and reliability.

Fig. 2 Accuracy verification of quantum mechanical calculations by comparing with experimental/calculation values. Reproduction of 55 (a) bond length and (b) bond angle values for 12 synthesized nitrogen–oxygen compounds, as well as 56 values for 11 compounds by prior reported calculations. Concrete bond length and angle data are provided in the ESI.†

3.2 Customized machine learning model: high accuracy and robustness

Among the six popular universal models evaluated, FM, SVM, and LR models exhibit the best performance in the binary classification of nitrogen–oxygen compound stability, as depicted by the Accuracy, F1, and Recall metrics in Fig. 3(a and b). The FM achieved the highest scores: Accuracy = 94.8, F1 = 96.2, and Recall = 98.7, with the best average scores being 93.5, 95.5, and 97.9, respectively. In comparison, LR model shows the highest scores: Accuracy = 94.3, F1 = 95.7, and Recall = 97.0, with corresponding average scores of 92.9, 94.6, and 96.0, respectively. The SVM model shows a slight decrease in the highest F1 scores by ∼0.7 compared to LR model, while remaining 97.1% in the highest Recall scores. Simultaneously, the average Accuracy and F1 scores for this model are slightly lower by 0.8 and 1.7 points, compared to the FM model, respectively, whereas its Recall scores show an improvement of 0.4. Overall, the FM, SVM, and LR models outperform the other three evaluated models-NB, GBDT, and RF-with their highest and average Accuracy scores being ∼1.2 and ∼5.0 points higher, respectively.

Fig. 3 Evaluation of accuracy and robustness of our customized machine learning model against various data splits. (a) The highest and (b) average performance scores of our customized model compared to those of six universal models in terms of accuracy. The performance metrics include Accuracy, F1, and Recall, represented by red, blue, and orange bars, respectively. (c) Statistics of Accuracy scores. Δ_AS denotes the Accuracy fluctuation, measured by the difference between the highest and lowest scores. The red and blue-filled frames delineate the central 50% of the data. The upper and lower curves connect the highest and lowest scores of the individual models. The black open squares represent average values, the horizontal solid lines indicate median values, and the inserted black lines are the error bars. The pie chart in the lower left corner depicts the composition of the customized model.

Based on these findings, we created a customized model by leveraging the advantage of FM, SVM, and LR models in handling small datasets and outliers. To determine the optimal mixing weight, we started with an FM : SVM : LR ratio of 5 : 0 : 5 and gradually adjusted the proportions. Inclusion of the SVM model significantly improved the metrics, and after extensive testing, the optimal ratio was found to be FM : SVM : LR = 4 : 3 : 3. This combination outperformed any individual universal model, achieving the best metrics: Accuracy = 94.9, F1 = 96.3, and Recall = 98.7, with corresponding average scores of 93.7, 95.6, and 98.0, respectively, as shown in Fig. 3(a and b).

The robustness of our customized model was evaluated by its Accuracy scores across different data splits, compared to those of six universal models, as shown in Fig. 3(c). Each model underwent 500 training cycles, with scores averaged over five-fold cross-validation. The index Δ_AS, representing the accuracy fluctuation between the highest and lowest scores, was used to characterize robustness. As depicted in Fig. 3(c), the 100 data points of our customized model fluctuated closely around the average value, with the lowest Δ_AS value of 2.81. This fluctuation is 34.5–280.0% smaller than those of the universal models NB, GBDT, and RF, and also surpassed its individual components, out-performing FM by 13.8%, SVM by 11.0%, and LR by 3.91%. Analysis of the error bar values, shown in black, reveals that the customized model has the smallest error bar, measuring only 0.0624, which reflects minimal fluctuation in performance across different data splits. More details are provided in Fig. S1 in the ESI.†

From above, our customized model surpasses universal models in both accuracy and robustness, achieving stable and high performance on small datasets of nitrogen–oxygen compounds. The model features a user-friendly graphical interface, making it easily accessible to researchers, as shown in Fig. S3 (ESI†).

3.3 Regulatory mechanism of stability in terms of energetics, electronic, and bonding

3.3.1 Thermodynamic stability: optimal ranges of energetics and electronic features. The 7820 energetics, electronic, and bonding parameters of the 168 nitrogen–oxygen compounds, calculated using high-accuracy quantum mechanical methods, were classified into two categories based on their thermodynamic stability. The Pearson correlation coefficient was calculated prior to machine learning. Fig. 4(a) shows that Pearson's correlation coefficients for all selected features are below 0.6, indicating weak correlations among the chosen features and confirming the validity of our feature selection.

Fig. 4 Regulatory mechanism of thermodynamic stability of nitrogen–oxygen compound in terms of energetic, electronic, and bonding. (a) Pearson's correlation coefficient and feature importance analysis for the four features. (b) Stability dependence on energy gap and dipole moment. The frames and curves represent statistics of the quantum mechanical calculation data, with the former showing the range of the central 50% of the data and the latter depicting the normal distribution of the full data. The scattered points represent the calculated data for 168 manually designed structures. Horizontal dashed lines mark each average value, corresponding to the peak of the normal distribution.

Feature importance analysis was performed after machine learning. Through feature importance analysis, we identified the most influential factors in regulating nitrogen–oxygen compound stability: energy gap, dipole moment, Laplacian bond order, and proportion of O-connected bond, in decreasing order of importance. The energy gap and dipole moment emerged as the dominant factors, with 74% and 16% importance percentages, respectively. Fig. 4(b) illustrates the distributions of these two factors. For both thermodynamically stable and unstable categories, the scatter plots show distinct aggregations around different average values, which correspond to the peak positions of their respective normal distributions. The significant differences in the average values between these two categories allow us to determine the optimal energy gap and dipole moment ranges for stable nitrogen–oxygen compounds. Based on the first and third quartiles of the stable category, the optimal energy gap range is > 13.46 eV, with higher energy gaps implying more difficult electron transitions and greater stability. The pivotal role of the energy gap in determining material stability has been extensively supported in previous studies.⁶² Similarly, the optimal dipole moment is < 1.54 Debye, as lower dipole moments indicate more uniform electron distribution and greater stability in nitrogen–oxygen compounds.

As the bonding features depend on local environmental factors, they act to fine-tune the thermodynamic stability of nitrogen–oxygen compounds. As illustrated in Fig. S2 (ESI†), these bonding features present small fluctuations and do not exhibit clear threshold values like energy gap and dipole moment.

3.3.2 Mechanical stability: fine-tuning of firmness via bonding features. The mechanical stability was characterized by bonding features such as Laplacian bond order and the proportion of O-connected bonds, which play an important role in characterizing bonding connections from the perspective of firmness. Fig. 5(a) illustrates three pairs of nitrogen–oxygen compounds with identical chemical formulas or geometrical configurations-e16vs.e15, b30vs.e14, and d2vs.d7-yet, the contrasting parties exhibit completely different conformations, bonding features, and ultimately, distinct levels of firmness. From the conformation perspective, e16 and e15 both share the formula N₆O₃ and possess a hexagonal nitrogen backbone, but the outermost oxygen in e15 is rotated by 30°. This slight rotation causes a significant alternation in bonding features, leading e15 (1.40) to have four times higher firmness than e16 (0.27). The essential factor behind this enhancement is the presence of electron delocalization in e15, visualized through orbital iso-surfaces (0.03) for the lowest π orbital and the nonbonding σ_LP orbital contributed by oxygen LP electrons. Comparing b30 and e14, both share the formula N₄O₆, yet they have drastically different dimensions of geometrical conformations. b30 adopts a nearly wave-like conformation, with a minimum Laplacian bond order of 0.07, while e14 has a highly symmetric adamantane-like cage structure,⁶³ yielding a much higher bond order of 0.44. Notably, such enhancement is attributed to the hyperstatic constraints in a cage-like structure, which possesses more constraints than necessary to maintain stability, making the firmness of e14 five times greater than that of b30. Lastly, d2 and d7 share the same pentagonal geometrical configuration, but their differing oxygen proportions lead to subtle changes in electron delocalization, resulting in d7 being three times firmer than d2. This demonstrates that introducing electron delocalization and hyperstatic constraint can efficiently enhance the firmness of nitrogen–oxygen compounds.

Fig. 5 Fine-tuning mechanical firmness of nitrogen–oxygen compound via bonding properties. (a) Illustration of three pairs of compounds with identical chemical formulas or geometrical configuration, but different firmness, and molecular orbitals: e16vs.e15, b30vs.e14, and d2vs.d7. The weakest bonds are highlighted in luminous magenta. Electron delocalization is visualized through orbital iso-surfaces (0.03) for the lowest π orbital and nonbonding σ_LP orbital contributed by oxygen lone pair (LP) electrons. (b) Established physical models of firmness and (c) size by fitting into the data from the 120 stable nitrogen–oxygen compounds. Dashed curves are the established models. Scatters represent the data of the 120 stable compounds. Red circles represent compounds with delocalized oxygen LPs, and blue squares represent compounds with all oxygen LPs localized.

Applying the above-revealed mechanism to all 120 designed stable compounds established a physical model outlining firmness, as depicted in Fig. 5(b). Oxygen LP delocalization plays a central role in regulating firmness, with additional factors including the number of oxygen atoms, their proportion, bond order, and the overall dipole moment. These factors govern the minimum Laplacian bond order following a logarithm growth trend. A physical model for nitrogen–oxygen compound size was also established, as shown in Fig. 5(c). Compound size is defined as the radius from the molecule's centroid to the outermost oxygen atoms. LP delocalization remains the dominant regulatory factor, alongside the outermost oxygen atoms' number, proportion, and bond order. These factors regulate compound size according to a quadratic growth trend. The fitting parameters of these models are listed in Tables S2 and S3 in the ESI.† It's noteworthy that a larger compound size correlates with reduced firmness.

Notably, compounds a1–a12 exhibit exceptional firmness and smallest sizes, reasonably interpreting their successful synthesis and widespread application. Specifically, the experimentally validated compounds N₂O (a5), and N₂O₅ (a12) exhibit bond firmness values of 1.49, and 0.23, respectively-all higher than the 0.07 of N₂O₄ (a8), which is widely used as an oxidizer in rocket propellants. Consistently, their experimentally measured decomposition temperatures are 600 °C, and 45 °C, respectively-all exceeding the 30 °C of N₂O₄ (a8). This consistency supports the alignment of our predicted mechanical stability with experimental trends. In particular, compound e15 exhibits a superior firmness of 1.40, which is comparable to the 1.49 firmness of the known stable compound N₂O (a5). As shown in Table 1, compounds c14, d7, d18, and e14 also demonstrate notable firmness values of 0.54, 0.52, 0.51, and 0.45, respectively. The firmness of tetra-nitrogen dianion (c14) and pentazole anion (d18) has been predicted in our prior studies, and the way to stabilize them via acidic entrapment was thoroughly analyzed.^64,65

Table 1 Mechanical stability of these designed nitrogen–oxygen compounds compared to the reported synthesized references, along with their decomposition temperature (T_d)

This study designed N–O compd	Mechanical stability	Reported synthesized reference
		Compounds	T d (°C)	Firmness
e15	1.40	a5 (N2O)	600.0 (ref. 66)	1.49
c14	0.54	a12 (N2O5)	45.0 (ref. 67)	0.23
d7	0.52
d18	0.51
e14	0.45
61 novel compd	>0.23
106 novel compd	> 0.07	a8 (N2O4)	30.0 (ref. 68)	0.07

---

Among the designed nitrogen–oxygen compounds, 106 exhibit greater firmness than N₂O₄ (a8) and 61 are firmer than N₂O₅ (a12), as shown in Fig. 5. These results suggest that the designed compounds, in particular those with small sizes, have the potential to withstand external heat or stress and maintain their structural integrity under varying conditions, indicating their suitability for practical applications.

3.4 Potential use as oxidizer in propellent: high energy density of typical stable oxidizers

From a practical usability perspective, we selected the three stable nitrogen–oxygen compounds e15, d7, and e14, to assess their potential as oxidizers in propellent. The energy density of the stable compounds was first evaluated through their heats of formation, with higher positive values indicating higher energy storage capability. As shown in Table 2, the enthalpies of formation were calculated to be positive, which are 664.7 kJ·mol⁻¹ for e15, 669.9 kJ·mol⁻¹ for e14, and 336.8 kJ·mol⁻¹ for d7, respectively, all significantly higher than those of conventional oxidizers.

Table 2 Enthalpy of formation (ΔH_f,solid) and specific impulse of three typical stable nitrogen–oxygen compounds and referenced oxidizers studied in this work

Oxidizers	ΔHf,solid (kJ mol^−1)	I sp in formula (1) (s)	I sp in formula (2) (s)
e15	664.7	295.2	302.4
e14	669.9	293.6	305.4
d7	336.8	262.6	270.3
AP	−295.8	263.9	261.4
RDX	70.3 (ref. 69)	257.2	251.5
HMX	104.8 (ref. 69)	256.4	249.4
CL-20	397.8	255.9	249.7

---

Next, we evaluated the energy density of the stable compounds using I_sp, which characterizes their capability of generating external work and launch distance of oxidizers when forming into propellant formulas. We considered the two most commonly used mixtures: formula (1), comprising 71% oxidizer by weight, 20% Al fuel, and 9% hydroxyl terminated polybutadiene (HTPB); and formula (2), comprising 73% oxidizer, 13% Al fuel, and 14% HTPB. Fig. 6 demonstrates that the specific impulses of e15, e14, and d7 generally surpass those of the current top-performing AP, RDX, HMX, and CL-20 oxidizers, in both formulations. The specific pulses of the three oxidizers range from 270.3 to 305.4 s, exceeding the best-performing CL-20 (249.6 s) by over 8.3%. Notably, e15 and e14 achieve specific impulses of 302.4 and 305.04 s, representing increases of 15.7% and 16.8%, respectively, compared to the most widely used oxidizer AP.

Fig. 6 Energy density level of three stable nitrogen–oxygen compounds (e15, d7, and e14) compared to current top-performing oxidizers (AP, RDX, HMX, and CL-20). The energy density level is characterized by specific impulse when formulated into formula (1) (71 wt% oxidizer, 20 wt% Al, and 9 wt% HTPB) and formula (2) (73 wt% oxidizer, 13 wt% Al, and 14 wt% HTPB).

Finally, by using three typical stable nitrogen–oxygen compounds-e15, e14, and d7-as examples, we demonstrate that our designed novel compounds exhibit exceptional firmness, allowing them to withstand ambient and extreme conditions. Moreover, their high specific impulses surpass those of all currently used oxidizers, predominantly forming into environmentally benign N₂, CO₂, and water products upon reaction with hydrocarbon fuels. This combination of high stability and high-energy density underscores their strong potential for application in rocket propellants, warranting the exploration of feasible synthetic routes despite the associated challenges. Notably, handling nitrogen–oxygen compounds could be risky, and researchers must undergo professional training and implement necessary protective measures to ensure safety.

4. Conclusions

To conclude, we developed a sequential construction strategy to design 168 nitrogen–oxygen compounds with distinct structural innovation and created a customized machine-learning model that demonstrated superior accuracy and robustness in predicting compound stability. The input feature descriptors were derived from high-throughput quantum mechanical calculations of 7820 structural and property parameters with the highest possible accuracy, and the reliability of the employed calculation method was thoroughly validated through comparison with experimentally characterized structural parameters. Our customized machine-learning model achieved the highest accuracy score of 94.9, outperforming the universal models. More importantly, the fluctuation in accuracy scores exhibits 13.8% greater robustness across various data splits, achieving stable and high performance on small datasets of nitrogen–oxygen compounds.

Data-driven analysis revealed the energy and electron-related characteristics as key factors in regulating the thermodynamic stability of nitrogen–oxygen compounds, while physics-driven insights uncovered that electron delocalization and hyperstatic constraints fine-tune their mechanical firmness. In this context, more constraints enable the compound to maintain stability under external forces. Among the designed nitrogen–oxygen compounds, 106 are expected to be more stable than synthesized compound N₂O₄, and 61 are more stable than N₂O₅, allowing them to potentially withstand ambient and extreme storage and application conditions. Simultaneously, the novel compounds are predicted to possess high specific impulses up to 305.4 s, surpassing all currently used oxidizers. This high stability and high-energy density combination highlights their potential for rocket propulsion.

Data availability

The data supporting this article have been included in the ESI.†

Conflicts of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

This work is supported by the National Natural Science Foundation of China (Grant No. 12222204 and 12472361) and the Foundation of State Key Laboratory of Explosion Science and Safety Protection (Grant No. QKKT24-01).

References

1. J. M. Wallace and P. V. Hobbs, Atmospheric Science, an Introductory Survey, Elsevier, Amsterdam, 2nd edn, 2006.
2. A. Snis and I. Panas, Mol. Phys., 1997, 91, 951–962.
3. Z. Sun, Y. D. Liu, C. L. Lv and R. G. Zhong, J. Mol. Struct.:THEOCHEM, 2009, 908, 107–113.
4. L. P. Cheng and W. Q. Cao, J. Mol. Model., 2007, 13, 1073–1080.
5. H. M. Badawi, W. Forner and K. S. Al-Ghamdi, Spectrochim. Acta, Part A, 2004, 60, 941–949.
6. L. P. Olson, K. T. Kuwata, M. D. Bartberger and K. N. Houk, J. Am. Chem. Soc., 2002, 124, 9469–9475.
7. R. Gilardi, J. Flippen-Anderson, C. George and R. J. Butcher, J. Am. Chem. Soc., 1997, 119, 9411–9416.
8. A. Snis and I. Panas, Chem. Phys., 1997, 221, 1–10.
9. J. K. Park and H. Sun, Chem. Phys., 2001, 263, 61–68.
10. T. M. Klapötke, A. Schulz and I. C. Tornieporth-Oetting, Chem. Ber., 1994, 127, 2181–2185.
11. R. D. Harcourt, J. Mol. Struct.:THEOCHEM, 1995, 342, 51–57.
12. E. D. Glendening and A. M. Halpern, J. Chem. Phys., 2007, 127.
13. A. Schulz, I. C. Tornieporth-Oetting and T. M. Klapötke, Angew Chem. Int. Ed. Engl., 1993, 32, 1610–1612.
14. Y. C. Hou, A. Janczuk and P. G. Wang, Curr. Pharm. Des., 1999, 5, 417–441.
15. E. Culotta and D. E. Koshland, Science, 1992, 258, 1862–1865.
16. W. J. You, H. J. Moon, S. P. Jang and J. K. Kim, Int. J. Heat Mass Transfer, 2013, 66, 101–110.
17. F. Zhang, H.-y. Chen, J.-c. Feng and D. Zheng, Fuel, 2021, 288, 119688.
18. R. Ahlrichs and F. Keil, J. Am. Chem. Soc., 1974, 96, 7615–7620.
19. Y.-i. Izato, K. Shiota and A. Miyake, Combust. Flame, 2021, 229, 111389.
20. W. R. Angus, R. W. Jones and G. O. Phillips, Nature, 1949, 164, 433–434.
21. L. E. Edwards, B. M. Kariuki, M. Didsbury, C. D. Jones and T. Wirth, Chem. Commun., 2024, 60, 5920–5923.
22. J. Lang and M. A. Bohn, Propellants, Explos., Pyrotech., 2021, 46, 233–244.
23. N. Wiberg, A. F. Holleman and E. Wiberg, Holleman-Wiberg's Inorganic Chemistry, Elsevier Science, 2001.
24. L. Zhu, Z. Wang, Y. Wang, G. Zou, H. K. Mao and Y. Ma, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 751–753.
25. T. M. Klapötke, Chemistry of High-Energy Materials, De Gruyter, 2015.
26. J. Zhang, J. Zhang, J. Singh, W. Wu, R. J. Staples, J. Zhang and J. n. M. Shreeve, J. Mater. Chem. A, 2024, 12, 5918–5923.
27. M. Rahm, S. V. Dvinskikh, I. Furó and T. Brinck, Angew. Chem., Int. Ed., 2011, 50, 1145–1148.
28. E. L. Yuan, J. I. Slaughter, W. B. Koerner and F. Daniels, J. Phys. Chem., 1959, 63, 952–956.
29. N. Wang, Y. Zhang, W. Wang, Z. Ye, H. Chen, G. Hu and D. Ouyang, Adv. Drug Delivery Rev., 2023, 196, 114772.
30. H. Wang, J. Feng, Z. Dong, L. Jin, M. Li, J. Yuan and Y. Li, npj Comput. Mater., 2023, 9, 200.
31. X. Guo, Z. Wang, J.-H. Yang and X.-G. Gong, J. Mater. Chem. A, 2024, 12, 10123–10136.
32. M. Han, J. F. Joung, M. Jeong, D. H. Choi and S. Park, ACS Cent. Sci., 2025, 11, 219–227.
33. G. Huang, F. Huang and W. Dong, Chem. Eng. J., 2024, 492, 152294.
34. N. Singh, K. Samanta, S. K. Maharana, K. Pal, S. Tretiak, A. Talapatra and D. Ghosh, J. Mater. Chem. A, 2025, 13, 9192.
35. S. Song, F. Chen, Y. Wang, K. Wang, M. Yan and Q. Zhang, J. Mater. Chem. A, 2021, 9, 21723–21731.
36. Y. Wang, Y. Liu, S. Song, Z. Yang, X. Qi, K. Wang, Y. Liu, Q. Zhang and Y. Tian, Nat. Commun., 2018, 9, 2444.
37. Y. Wen, L. Wen, B. Tan, J. Dou, M. Xu, Y. Liu, B. Wang and N. Liu, J. Mater. Chem. A, 2024, 12, 9427.
38. R. Tsyshevsky, P. Pagoria, A. S. Smirnov and M. M. Kuklja, J. Phys. Chem. C, 2017, 121, 23865–23874.
39. J. Liu, S. Zhao, B. Duan, X. He, C. Yang, X. Pu, X. Zhang, Y. Xiao, F. Nie, W. Qian, G. Li and C. Zhang, J. Mater. Chem. A, 2023, 11, 25031–25044.
40. H. Li, L. Zhang, N. Petrutik, K. Wang, Q. Ma, D. Shem-Tov, F. Zhao and M. Gozin, ACS Cent. Sci., 2020, 6, 54–75.
41. L. Wen, Y. Wang and Y. Liu, Acc. Mater. Res., 2025, 6, 64–76.
42. W.-L. Yuan, L. Zhang, G.-H. Tao, S.-L. Wang, Y. Wang, Q.-H. Zhu, G.-H. Zhang, Z. Zhang, Y. Xue, S. Qin, L. He and J. n. M. Shreeve, Sci. Adv., 2020, 6, eabb1899.
43. X. Huang, C. Li, K. Tan, Y. Wen, F. Guo, M. Li, Y. Huang, C. Q. Sun, M. Gozin and L. Zhang, iScience, 2021, 24, 102240.
44. Z. Zhang, C. Chen, Y. Cao, L. Wen, X. He and Y. Liu, Thermochim. Acta, 2024, 735, 179717.
45. J. Wu, S. Song, X. Tian, Y. Wang and X. Qi, Energ. Mater. Front., 2023, 4, 254–261.
46. M. J. Frisch, G. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci and G. A. Petersson, Gaussian 09W, revision A.02, 2009.
47. T. Lu and F. Chen, J. Comput. Chem., 2012, 33, 580–592.
48. G. Alibaba, https://github.com/alibaba/Alink.
49. A. Hassanpour, S. Yasar, A. Ebadi, S. Ebrahimiasl and S. Ahmadi, J. Mol. Model., 2021, 27, 124.
50. M. Arhangelskis, A. D. Katsenis, N. Novendra, Z. Akimbekov, D. Gandrath, J. M. Marrett, G. Ayoub, A. J. Morris, O. K. Farha, T. Friščić and A. Navrotsky, Chem. Mater., 2019, 31, 3777–3783.
51. X. Hao, J. Liu, H. Luo, Y. Han, W. Hu, J. Liu, J. Li and X. He, Crystals, 2019, 9, 256.
52. R. Fournier and A. Mohareb, J. Chem. Phys., 2016, 144.
53. D. Han, T. Zhang, M. Huang, D. Sun, M.-H. Du and S. Chen, APL Mater., 2018, 6, 084902.
54. G. De Bo, Chem. Sci., 2018, 9, 15–21.
55. Z. Fu, H. Zhang, C. Si, D. Legut, T. C. Germann, Q. Zhang, S. Du, J. S. Francisco and R. Zhang, J. Phys. Chem. C, 2018, 122, 4710–4722.
56. T. Lu and F. Chen, J. Phys. Chem. A, 2013, 117, 3100–3108.
57. T. Lu and F. Chen, J. Comput. Chem., 2012, 33, 580–592.
58. S. M. J. Rogge, M. Waroquier and V. Van Speybroeck, Acc. Chem. Res., 2018, 51, 138–148.
59. G. Pavlovets and V. Tsutsuran, Physicochemical Properties of Powders and Propellants, Russian Ministry of Defense Publishing House, Moscow, Moscow, 2009.
60. I. L. Dalinger, K. Y. Suponitsky, T. K. Shkineva, D. B. Lempert and A. B. Sheremetev, J. Mater. Chem. A, 2018, 6, 14780–14786.
61. P. Politzer, J. S. Murray and C. Eur, J. Energ. Mater., 2011, 8, 209–220.
62. N. R. Badders, C. Wei, A. A. Aldeeb, W. J. Rogers and M. S. Mannan, J. Energ. Mater., 2006, 24, 17–33.
63. M. Bremer, P. von Ragué Schleyer, K. Schötz, M. Kausch and M. Schindler, Angew. Chem., Int. Ed., 1987, 26, 761–763.
64. L. Zhang, C. Yao, Y. Yu, S.-L. Jiang, C. Q. Sun and J. Chen, J. Phys. Chem. Lett., 2019, 10, 2378.
65. C. Jiang, L. Zhang, C. Sun, C. Zhang, C. Yang, J. Chen and B. Hu, Science, 2018, 359, eaas8953.
66. G. Löffler, V. J. Wargadalam, F. Winter and H. Hofbauer, Combust. Flame, 2000, 120, 427–438.
67. P. Connell and H. S. Johnston, Geophys. Res. Lett., 1979, 6, 553–556.
68. V. Ragaini, Sol. Energy, 1982, 29, 535–540.
69. Y. Tang, C. He and J. n. M. Shreeve, J. Mater. Chem. A, 2017, 5, 4314–4319.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5ta00267b

---