DeepEMs-25: a deep-learning potential to decipher kinetic tug-of-war dictating thermal stability in energetic materials

Abstract

Atomic-scale insight into decompositions in energetic materials (EMs) is essential for harnessing energy release, which remains elusive due to both instrumental and computational limitations. Herein, we developed DeepEMs-25, a deep-learning potential trained on diverse EMs towards accurate and efficient simulations. Applying DeepEMs‑25 to an isostructural ABX₃ molecular perovskites series, with A-site organic cations, B-site alkali or ammonium cations, and X-site perchlorate anions, we probe the effect of cation size on reactivity. Arrhenius analysis of 100-ps trajectories reveals that increasing B‑site ionic radius simultaneously decreases X–A collision’s activation energy (enhancing reaction rates) and decreases X–A collision’s pre‑exponential factor (reducing collision frequency), producing opposing kinetic effects. Such “kinetic tug‑of‑war” explains why an intermediate‑sized cation yields maximal thermal stability by optimally balancing reactivity and collision dissipation. A similarly sized reactive cation promotes additional hydrogen-transfer pathways causing accelerating decomposition. Our findings link atomistic kinetics to macroscopic stability, informing next-generation EMs design.

Introduction

Energetic materials (EMs) epitomize humanity’s quest to harness chemical energy at atomic scales^{1, 2, 3, 4, 5, 6, 7–8}. The development of advanced EMs is driving a paradigm shift from traditional organic synthesis toward modular, multi-ion assembly strategies. Energetic molecular perovskites in ABX₃ frameworks where organic A‑site (typically larger organic cations) and B-site (typically smaller inorganic cations) cooperate with X-site oxidizing anions (e.g., ClO₄^–) to build oxidizer-fuel integrated structures, have attracted much attention for their controllable energy release and thermal stability (Fig. 1a). However, experimental studies have found that the macroscopic performance changes of isomorphous series are non-monotonic. For instance, in the representative series, 1,4-diazabicyclo[2.2.2]octane-1,4-diium (H₂dabco²⁺) alkali/ammonium perchlorate series (DAPs), e.g., (H₂dabco)B(ClO₄)₃, (B = Na⁺, K⁺, Rb⁺, and NH₄⁺ for DAP-1, DAP-2, DAP-3, and DAP-4, respectively, Fig. 1b)^{9, 10, 11–12}, B-site ionic radii increase monotonically following the well-established periodic trend, Na⁺ < K⁺ < NH₄⁺ < Rb⁺, which would simply suggest a corresponding monotonic trend in macroscopic properties. However, the observed behavior significantly diverges from this anticipation that K⁺-based DAP-2 exhibits peak thermal stability (onset decomposition temperature T_d = 364 °C) and friction sensitivity (42 N), whereas NH₄⁺-based DAP-4 shows abrupt decomposition despite the similar ionic radii of NH₄⁺ and K⁺¹². These anomalies underscore the inadequacy of static microscopic descriptors and highlight the need for atomistic interrogation of multi-factor and dynamic coupling.

Fig. 1 Structures and possible reaction pathways in DAPs. [Images not available. See PDF.]

a The schematic template of molecular perovskite. b The unit cells of DAP-1, DAP-2, DAP-3, and DAP-4. All possible initial inter-ionic decomposition paths of energetic molecular perovskites with unreactive (c) and reactive (d) B-site ions.

Traditional methods face challenges in analyzing such rapid and complex solid-state reactions. First, experimental techniques such as differential scanning calorimetry (DSC), infrared spectroscope (IR), and thermogravimetric analysis (TGA), lack the spatiotemporal resolution needed to capture initial bond‑breaking events, yielding only bulk decomposition profiles^{13, 14, 15, 16–17}. Second, computational tools face complementary limitations. (i) density functional theory (DFT) is restricted to hundreds of atoms and is difficult to capture statistical decomposition events¹⁸; (ii) the empirical reaction force field (ReaxFF) has low parameter transfer ability and is difficult to cover emerging energetic molecular perovskite systems^19,20. Likewise, the state-of-the-art datasets for deep learning potentials, although excellent in organic gas-phase chemistry (e.g. Transition-1x²¹, ANI-1xBB²², and OMol25²³) or CHNO organic EMs²⁴, they fall behind the pace with the rapid emergence of new EMs and lack widespread solid-state reaction dataset support.

To address the challenges, we developed DeepEMs‑25, a deep-learning potential (DP) trained on a dataset covering decomposition trajectories for 20 EMs, extending beyond conventional CHNO organic EMs to include oxoanion salts and azide salts spanning 13 elements, where the completeness was ensured by active learning. Built within the DeePMD framework^{25, 26–27}, DeepEMs‑25 delivers high accuracy comparable to DFT while scaling linearly with system size, enabling large‑scale molecular dynamics of solid‑state reactions without biasing methods.

Applying DeepEMs-25, we investigated the four highly isostructural energetic molecular perovskites on initial decomposition behavior. By systematically enumerating all possible initial inter-ionic reaction pathways in the DAP series, we identified two distinct scenarios: (i) pathways dominated by X–A interactions when site B remains unreactive (Fig. 1c), and (ii) pathways involving active participation of site B, which can significantly alter the decomposition kinetics (Fig. 1d). By focusing on those dominant pathways, we revealed an antagonistic reactivity mechanism involving (i) X–A attraction and (ii) B-site steric effects, which dictates the decomposition behavior. This discovery not only explains the counterintuitive thermosensitive behavior of the DAP series, but also provides a new perspective for the study of solid-state chemical reaction mechanisms.

Results

Development and evaluation of DeepEMs-25

To enable reliable, large‑scale simulations of solid‑state reactions, we developed DeepEMs-25, a deep-learning interatomic potential^28,29 trained via DeePMD framework^26,27,30 which decomposes the total energy E into per‑atom contributions and recovers forces F and virials V by analytic differentiation, thereby ensuring strict energy conservation and stable long‑trajectory dynamics. The dataset explicitly covers 20 EMs (Supplementary Fig. S1), including organic explosives, oxoanion salts, and azide salts. Progressive active learning^{31, 32–33} sampled diverse configurations up to 4000 K and 120 GPa (Supplementary Table S1) and labeled with high-quality spin-polarized DFT labels guaranteed by rigorous convergence (Supplementary Figs. S2–S4), yielding a final dataset with over 87000 frames (Supplementary Figs. S5, S6). We visualized the model’s learned chemical representations using t‑distributed stochastic neighbor embedding (t‑SNE) on its latent atomic descriptors (Fig. 2a), observing well‑defined clusters that match established chemical families (Fig. 2b).

Fig. 2 The sketch-map visualization of the DeepEMs-25 datasets. [Images not available. See PDF.]

a The t-SNE visualizations of the learned descriptor of the DP model and b sketch crystal structures visualizations for 20 EMs covered by our DP model. See Supplementary Fig. S1 for further details.

Before formal MD running, we conducted extensive evaluations to validate its performance. A high-quality DP model requires (i) high accuracy in directly predicting DFT labels, (ii) high efficiency running inferences, (iii) conservation in potential energy surface, and (iv) reliable domain-specific downstream property calculations. First, the DeepEMs-25 model demonstrates ca. 100 times faster than DFT-MD and superior to semi-empirical PM6 (Fig. 3a). Second, direct inference on an independently constructed test set (see Methods for details) yielded root mean square errors (RMSE) of 26.6 meV/atom for energy, 375.7 meV/Å for force, and 21.7 meV/atom for the virial coefficient (Fig. 3b). Third, the conservativeness of the yielded DP models was confirmed by evaluating energy drifts and standard deviations in NVE-MD, thereby validating its MD simulation reliability (Supplementary Fig. S7). Furthermore, to quantify early decomposition events in crystals, we introduce an in‑lattice bond‑cleavage benchmark. Under periodic boundary conditions, we track the energy-distance curve of key bonds of N–H in ammonium, N–H, C–H, C–N in H₂dabco, and Cl–O in perchlorate, using structures generated by systematically stretching those bonds from 0.75 to 1.5 times their equilibrium lengths in 12 incremental steps and compare them between DeepEMs-25 and DFT. DeepEMs‑25 reproduces the DFT cleavage curves with high fidelity, establishing this metric as a robust indicator of bond‑breaking propensity within confined lattices (Fig. 3c).

Fig. 3 Evaluations on DeepEMs-25 model. [Images not available. See PDF.]

a Efficiency evaluations, where CPU tests were conducted on Intel Xeon Platinum 8358P CPU @ 2.60 GHz and GPU tests were conducted on NVIDIA RTX 4090, 24GB and NVIDIA V100, 32GB, and b direct inference on test set results compared to DFT, where the color mapping represents the kernel density estimation using Scott’s bandwidth rule and equal weighting. c The in-lattice bond cleavage benchmarks, with the inset illustrating the bond stretching schematic for a representative species within the lattice.

Collectively, those results demonstrate that DeepEMs‑25 combines exceptional computational efficiency, DFT‑level accuracy, and rigorous energy conservation, enabling high‑fidelity, large‑scale simulations for EMs.

Species and reaction analysis for DAPs decomposition

Having established DeepEMs-25, we next apply it to investigate the initial decompositions of DAPs. To elucidate the decomposition mechanisms of the DAP series, we conducted MD simulations in the temperature range of 1666–2500 K using 3 × 3 × 3 supercells of DAP-1 to DAP-4. The simulations were performed in the NVT ensemble after energy minimization and NPT equilibration at 300 K and 1 bar. To monitor the stability of the MD process on the fly, we calculated the average deviations for forces, which were consistently below 5% (Supplementary Fig. S8). This metric served as a lightweight surrogate for simulation robustness, given that reference DFT forces are impractical to compute for large supercells and long trajectories.

Accurate identification and tracking of chemical species during decomposition are critical for understanding the reaction mechanisms. To this end, two distinct methods were employed to track species evolution during decomposition. The first method is that for simple inorganic substances, it can be directly defined by the coordination number and cutoff radius outside the central atom (the criteria shown in Supplementary Table S2). For complex organic substances, it is necessary to further search the chemical species using a depth-first search algorithm, which is implemented by calling the OpenBabel library³⁴ in ReacNetGenerator by Zeng et al.³⁵ (the criteria shown in Supplementary Table S3) and recorded in SMILES encoding.

We define “the completion of initial reactions” as the reaction between site X and site A. That is, in this scenario, the reactions involving perchlorate anions (X) and H₂dabco²⁺ organic cations (A) are completed, so that the statistical number of remaining X and A in the system tends to 0 (Supplementary Figs. S9–S12). Through the above statistical method, we found that at 1666 K, perchlorate (X) and H₂dabco²⁺ (A) can be completely reacted within 100 ps. Furthermore, by tracking the species evolution of each atom in different frames, the reaction path can be constructed³⁵.

First, we investigate the dominant decomposition pathway of H₂dabco²⁺, focusing on their behavior within various DAP variants. Our statistical analysis reveals that the primary reactions involve two inter-ionic hydrogen transfer reactions and one intra-ionic ring-opening reaction: (i) deprotonation of N–H and hydrogen transfer between perchlorate and H₂dabco²⁺, (ii) bond homolysis of C–H and hydrogen transfer between perchlorate and H₂dabco²⁺, and (iii) cleavage of the C–N bond and ring opening. Consistently, across all four DAP variants examined, the hierarchy of these pathways is maintained, with N–H deprotonation emerging as the most prominent mechanism, followed by C–H homolysis, and finally C–N ring opening (Fig. 4a–d). At relatively low simulation temperatures (e.g., 1666 K), more than 60% of the H₂dabco²⁺ in all four materials undergoes N–H deprotonation, about 30% undergoes C–H homolysis, and C–N ring opening accounts for at most 5%. As the temperature increases, the distribution of products shifts significantly, with the reaction ratios for C–H cleavage and C–N bond cleavage increasing. For instance, at 2500 K, the proportion of C–N ring opening rises to ~20%. Nevertheless, the reactions remain predominantly governed by N–H and C–H cleavage.

Fig. 4 Temperature-dependent decomposition products evolution of H₂dabco²⁺. [Images not available. See PDF.]

Transient counts of key decomposition species from H₂dabco²⁺ at various temperatures in a DAP-1, b DAP-2, c DAP-3, and d DAP-4. Insets are percentages of cumulated species counts.

To gain deeper insights, we employed the protocol by Wu et al.³⁶ to calculate the reaction rates of the chemical reactions corresponding to these three pathways at different temperatures and performed Arrhenius fitting (Fig. 5a). The results are consistent with chemical intuition: the activation energy of the N–H bond is in the range of 55.69 ± 5.40 to 64.10 ± 8.17 kJ/mol (Fig. 5b), because N–H deprotonation is an acid-base reaction; the activation energy of the C–H bond cleavage is the highest, ranging from 109.64 ± 2.69 to 120.29 ± 6.50 kJ/mol (Fig. 5c); so the latter is sensitive to temperature and its proportion increases at high temperatures. The activation energy of the C–N bond cleavage is in the range from 68.64 ± 2.98 to 72.11 ± 2.82 kJ/mol (Fig. 5d). The significant temperature sensitivity of the C–N cleavage reaction implies a high activation entropy. Upon ring opening, the liberated chain segment gains conformational freedom, increasing the number of accessible spatial configurations. Simultaneously, the electronic structure transitions from a restricted closed-shell state to an open-shell diradical state. Both effects result in a large entropic contribution to the entropic barrier, rendering the reaction highly temperature-sensitive despite its relatively modest activation energy.

Fig. 5 Kinetic analysis on X–A reaction pathways. [Images not available. See PDF.]

a Schematic illustration of decomposition pathways of H₂dabco²⁺. The Arrhenius fitting of hydrogen transfers on b N–H and c C–H, where the shades denote the 95% confidence interval. d The relevant fitted activated energies and pre-exponential factors for different materials. Comparison of fitted activation energies and pre-exponential factors for DAPs for e hydrogen transfer involving N–H, f hydrogen transfer involving C–H, and g ring opening on C–N. Units for pre-exponential factors A and rate constants k are 1/s for C–N cleavage and cm³/(mol·s) for hydrogen transfer, omitted for clarity. The representative MD snapshots from DAP-4 at 2000 K of h deprotonation on N–H, i homolysis on C–H, and j C–N cleavage.

Next, we study the regularity of the kinetic parameters of different materials in each elementary reaction. Comparing DAP-1, DAP-2, and DAP-3, we find that with the increase of the radius of the cation at the B site, the activation energy and pre-exponential factor of the C–H homolysis and N–H deprotonation reactions decrease respectively; comparing DAP-2 and DAP-4, we find that the activation energy and pre-exponential factor are close (Fig. 5e–f). For ring opening on C–N, the difference is relatively not prominent (Fig. 5g). This means that the reaction kinetic parameters at the X and A sites are highly correlated with the properties of the B site, mainly the radius. However, these two results have opposite effects on the X–A reaction. First, the reduction of activation energy lowers the barrier for X–A reactions, while the decrease in the pre-exponential factor indicates fewer effective collisions. This antagonistic relationship between E_a and A has also been observed and investigated in thermally activated solid-state processes^37,38. In other words, although the monotonic B-site properties point to two monotonic changes in kinetic parameters, the effects of these two kinetic parameters on the apparent reactivity are opposite. This likely leads to the optimal balance of macroscopic thermal decomposition temperature and friction sensitivity.

For the C–N ring opening pathway, the relatively minor differences in the activation energy and pre-exponential factor across DAP-1, DAP-2, and DAP-3 (Fig. 5g) may stem from the primarily intramolecular nature of this reaction and the comparatively weaker influence exerted by the B-site cation. Unlike N–H deprotonation and C–H homolysis, which are more sensitive to the electronic and steric effects of the B-site cation, C–N ring opening is primarily an intramolecular process involving the cleavage of the C–N bond within the H₂dabco²⁺ ring. This reaction is less dependent on the interactions between site X (perchlorate) and site A (H₂dabco²⁺), and thus the properties of the B-site cation (e.g., radius, charge density) have a weaker influence on its kinetics. This observation aligns with established solid-state reactivity observations³⁹ suggesting that intramolecular processes are generally less sensitive to lattice perturbations than intermolecular reactions. Consequently, its kinetic parameters are less affected by B-site variations.

Complementary to these kinetic analyses, representative MD trajectory snapshots captured at 2000 K (Fig. 5h–j) visually depict the hydrogen transferring processes on the N–H and C–H bonds as well as the C–N ring opening events at the atomic level. These snapshots highlight the key reactive events, with specific particle identifiers marking the species involved, thereby reinforcing the proposed reaction pathways.

Influence of static properties of site B

To further analyze those trends, we examined the non-chemical kinetic effects of B-site radius. Since site B is non-reactive in DAP-1, DAP-2, and DAP-3, “ineffective collisions” between the X site and B site compete with the effective collisions of X–A. We defined the collision frequency (C.F.) of X–B as follows:

1

$$\mathrm{C}.\mathrm{F}.=\frac{C_{\mathrm{B},\mathrm{coll}}}{N_{\mathrm{X}}*T_{\mathrm{sim}}}$$ where C_B,coll represents the counts of collision events for B, defined as those events where the distance between the X and B sites is less than a threshold determined from the ionic radii (set to 1.2 times (r_X + r_B) as listed in Supplementary Table S4); N_X represents the counts of site X, placed in the denominator for normalization; T_sim represents the simulation time span for statistics. As shown in Fig. 6a, C.F. in different time intervals at 1666 K shows that in the early stage of the reaction (simulation time < 0.1 ps), the collision frequency exhibits a pronounced positive correlation with the ionic radius. The observed trend was DAP-3 > DAP-4 ~ DAP-2 > DAP-1, suggesting that larger B-site cations increase the likelihood of ineffective collisions. This, in turn, thereby reduces the probability of the more productive X–A collisions that drive the reaction forward.

Fig. 6 Interaction analysis of site X and site B. [Images not available. See PDF.]

a Collision frequency (C.F.); b Basin analysis of electron density (isovalue = 1.0 a.u.); c SAPT decomposition results.Influence of Reactivity of Site B.

From the perspective of activation energy, the effect of the B site on the X–A reaction appears to be predominantly mediated by the attractive interactions between the X and B sites. To investigate this phenomenon, we employed electron density basin analysis^{40, 41–42} to characterize the spatial region occupied by the B-site cations. The electron density isosurfaces (Fig. 6b, isovalue = 1.0 a.u.) show that the extent to which the electron density basins correlate with the ionic radii follows the order Na⁺ < K⁺ ~ NH₄⁺ < Rb⁺. When the geometric center separation between the X site and the B site is minimal, a larger atomic volume for the B-site cation leads to an enhanced interaction with the electron cloud associated with the X site. This interaction is predominantly governed by electrostatic attraction and inductive forces, with minimal contributions from dispersion effects, and is further modulated by mutual repulsion due to orbital exchange.

To further validate those observations, energy decomposition analysis was performed using symmetry-adapted perturbed theory (SAPT)^43,44. The results revealed that the three attractive forces (electrostatic, inductive, and dispersion) between X and B are relatively similar across the four materials, with variations of less than 6 kJ/mol. The most significant divergence arises from exchange repulsion, which follows the order Na⁺ < K⁺ < NH₄⁺ < Rb⁺. Consequently, the overall binding energy, when interacting with perchlorate, adheres to the trend Na⁺ > K⁺ > NH₄⁺ > Rb⁺ (Fig. 6c). This observation is in strong agreement with the experimentally determined activation energies, thereby substantiating the hypothesis that the attractive X–B interaction acts as an impediment to the progress of the effective X–A reaction.

The reactivity of the B-site cation plays a pivotal role in dictating the decomposition pathways and kinetics of energetic molecular perovskites. By comparing DAP-2 and DAP-4, we observe that although the ammonium ion has a radius similar to K⁺, its reactivity as a hydrogen carrier significantly modifies the reaction pathway. In DAP-4, the ammonia molecule, formed through the deprotonation of the B-site ammonium ion (Fig. 7a), actively participates in the decomposition process, leading to a dramatic reduction in the activation energy of key reactions. Specifically, the activation energy for C–H bond cleavage with ammonia in DAP-4 is reduced to 51.76 ± 8.66 kJ/mol (Fig. 7b), which is almost half of that with the participation of perchlorate. This substantial decrease underscores the catalytic role of ammonia in facilitating hydrogen transfer and lowering the energy barrier for bond dissociation (Fig. 7c).

Fig. 7 Hydrogen transfer enhanced by reactive B-site ammonium. [Images not available. See PDF.]

a Schematic illustration of hydrogen transfer pathway of H₂dabco²⁺ interacted with ammonia of DAP-4. b the Arrhenius fitting of hydrogen transfer on C–H, where the shades denote the 95% confidence interval. Units for pre-exponential factors A and rate constant k are both cm³/(mol·s) for hydrogen transfer, omitted for clarity. c The representative MD snapshots for DAP-4 at 2000 K.

Further statistical analysis of ammonium ions and ammonia dynamics (Supplementary Fig. S12) provides additional insights into the reaction mechanism. During the early stages of decomposition, the concentration of ammonium ions decreases, while ammonia initially increases. This trend reverses in the later stages, with ammonium ions increasing and ammonia exhibiting a fluctuating pattern. Notably, the point at which ammonium ion concentration reaches its minimum coincides precisely with the peak in H₂dabco²⁺ C–H homolysis products. This correlation suggests that H₂dabco²⁺ acts as a hydrogen donor in the ammonia-ammonium equilibrium, consuming ammonia and driving the equilibrium toward ammonium generation.

This dual role of H₂dabco²⁺, both as a reactant and a hydrogen donor, underscores the complexity of the decomposition process in DAP-4. The ability of ammonia to participate in hydrogen transfer reactions not only accelerates bond cleavage but also introduces a dynamic equilibrium that influences the overall reaction kinetics. These findings further demonstrate that the reactivity of the B-site cation is not merely an intrinsic static property; rather, it is a dynamic factor that can profoundly reshape both the decomposition pathway and the associated EMs.

Discussion

In this work, we employed large‑scale, unbiased molecular dynamics powered by our DeepEMs‑25 deep potential to uncover a collision‑controlled reactivity mechanism in a family of molecular perovskites, i.e., (H₂dabco)B(ClO₄)₃ (DAP-1 to DAP-4), with an unprecedented large-scale MD simulation. Our findings reveal a broadly applicable antagonistic kinetic principle that governs the decomposition behavior of these materials, bridging the gap between microscopic dynamics and macroscopic properties. Importantly, leveraging a high-quality deep potential model enables the resolution of rare decomposition pathways, some accounting for only a few percent of events, which would be statistically inaccessible in small simulation cells due to insufficient sampling in conventional small-scale simulations (e.g. the ring-opening pathway observed below 2000 K). Such limitation for conventional methodology is exemplified by the DAP-4 unit cell, which comprises only eight (ABX₃) formulas; such a limited structural size underscores that the statistically significant trends and emergent phenomena observed here are accessible only via large-scale simulations. This observation resonates with Anderson’s “more is different” concept⁴⁵, highlighting that collective behavior at extended scales (the coupling of chemical reaction pathways in this case) gives rise to properties that cannot be simply extrapolated from isolated or minimal configurations.

As summarized in Table 1, the key insight lies in the dual role of the B-site ionic radius: increasing ionic radius (Na⁺ < K⁺ < Rb⁺) simultaneously reduces the activation energy E_a, which enhances reactivity and decreases the pre-exponential factor A, which suppresses collision probabilities. This antagonistic coupling explains the peak thermal stability and insensitivity of K⁺-based DAP-2, which optimally balances these opposing effects. In contrast, the ammonium-based DAP-4, despite its ionic radius parity with K⁺, exhibits distinct behavior due to the reactive nature of NH₄⁺, which facilitates hydrogen transfer and promotes C–H bond rupture, leading to abrupt decomposition.

Table 1. Summary of identified microscopic factors and their relationship with macroscopic decomposition properties

Material	B-site cation	B-site radius (Å)62,63	ln[A(X–A)]a,b	Ea(X–A) (kJ/mol)b	Ebinding(X–B) (kJ/mol)	Td (°C)12
DAP-1	Na+	1.02	36.31 ± 0.39	120.20 ± 5.20	–401.44	344
DAP-2	K+	1.38	35.92 ± 0.13	113.73 ± 2.15	–391.04	364
DAP-4	NH4+	1.48	35.74 ± 0. 20	112.73 ± 2.26	–390.02	358
DAP-3	Rb+	1.52	35.72 ± 0. 16	106.84 ± 2.69	–379.66	352

^aPre-exponential factors are in cm³/(mol·s) for H transfers (the same as rate constant k, omitted for clarity).

^bDisplayed values refer to H transfer from C on H₂dabco²⁺ to ClO₄^–.

By linking static descriptors (ionic radius, complexion energies) to dynamic parameters (activation barriers and collision frequencies), our results not only resolve the counterintuitive trends in the DAP series but also provide a cornerstone for the rational control of energy release especially for EMs with tailored thermal stability and reactivity.

Looking forward, this collision‑controlled paradigm could be extended to more complex solid‑state systems, offering a path toward predictive control of reaction kinetics. Integrating deep potentials with emerging solid‑state reaction theories may ultimately broaden our ability to engineer advanced functional materials with unprecedented precision.

Methods

Datasets preparation

The crystal structures are collected from the Cambridge Crystal Data Center (CCDC) or appendixes in documents (Supplementary Fig. S1).

The QuickStep module of the CP2K package⁴⁶ is used for DFT calculations in dataset labeling. Since the quality of DFT data is crucial, our research started with a series of convergence tests. Considering the balance of computational cost, we utilized proper supercell (see Supplementary Fig. S1 for details) to allow gamma point-only calculations with orbital transformation (OT)^47,48 approach for DFT calculations. So the convergence tests were mainly conducted for cutoff energies (Supplementary Figs. S2–S4) and finalized the parameters to be 800 Ry with relative cutoff of 60 Ry. During DFT calculations, unrestricted shell was considered to appropriately describe generated radicals with explicit spins. All atoms are described using the MOLOPT basis set in combination with norm-conserving Goedecker-Teter-Hutter (GTH) pseudopotentials. The Perdew-Burke-Ernzerhof (PBE)⁴⁹ exchange-correlation functional is employed throughout and the electron density is written over an auxiliary plane-wave basis set with appropriate energy cut-off. The D3 dispersion correction was employed to appropriately correct the inter-molecular non-covalent interactions^50,51. The self-consistent field (SCF) convergence of both the outer and inner loops is achieved with an accuracy of 10⁻⁵ Hartree.

The initial training dataset was generated from short DFT-MD simulations (0.1 fs timestep, 500 fs total) for each material at 300 K and 2000 K (NPT-F ensemble, 1 bar), and 3500 K (NVT ensemble). Data including structure, energy, force, and virial were recorded and down-sampled to avoid redundancy (1% of 300 K trajectories, 5% of others).

To evaluate model performance, two types of test sets were constructed. (i) The independent test set comprised new DFT-MD trajectories for 26 structures (including isomers), each run at 3500 K (NVT, 5000 steps), with 50 frames randomly down-sampled per trajectory (total 1300 frames). These trajectories were generated using distinct initial velocities from those used in training, and were excluded from both training and active learning. (ii) For random split evaluation, we partitioned the full active-learning dataset into 95/5 train/test splits, repeated across five trials with different seeds. Each model was retrained from scratch and evaluated using a consistent training schedule.

Active learning

In each iteration of the active learning, four different DP models were trained with different random seeds for initial weights of the neural networks. After this, a short MD simulation was performed based on the first DP model, and the potential energies and atomic forces were also predicted by the other three models. Because atomic forces in thermal decomposition span a wide dynamic range, the relative model deviation was employed instead of the absolute for selecting structures during the simulation. For each atom in a given snapshot, the relative model deviation σ_i of force is calculated as:

2

$${\sigma }_i=\frac{\epsilon }{⟨f_i⟩+d}$$ where ⟨f_i⟩ denotes the average forces on the atom $i$ of the four DP models, and ε denotes the standard deviation of the atomic force. d is a constant used to make sure that the atom with a small ε will also have a small relative model deviation, and it was set to 1.0 eV/Å in this study. To obtain a frame-level uncertainty indicator, the maximum deviation σ = max(σ_i) over all atoms within the frame was used. If σ is less than the lower trust level, the snapshot will be labeled as “accurate”, which means that it has been accurately learned by the neural network; while if the σ of a given snapshot is within a certain range (lower trust level < σ < higher trust level), this snapshot will be labeled as “candidate”; else the snapshot was labeled as “fail” due to the potential crashed simulation. In each iteration, we randomly selected at most 250 frames of “candidate” structures. The potential energy, atomic forces, and virials of those selected frames were calculated by DFT and added back to the dataset. In each iteration, we appropriately increased the length of the simulation as listed in Supplementary Table S1.

Model architecture and training

To develop DeepEMs-25, we employed the DeepPot-Attention architecture without the attention layer (noted as DPA1-L0)³⁰, implemented in the DeePMD-kit code package (version 2.2.9)^25,52,53.

Following the standard DPA1 architecture³⁰, the total potential energy E = ∑_ie_i, where e_i is the energy of atom i, which depends on its local environment defined by the relative positions of neighboring atoms within a cutoff radius R_c of 6.0 Å and smooth cutoff radius R_cs of 0.5 Å. Each atomic environment is first transformed into generalized coordinates ${\tilde{R}}_i$ which encode normalized distance vectors with radial weighting. Simultaneously, the atomic type a_i is mapped into a length-fixed vector $T_i\in {\mathbb{R}}^8$ . The matrices ${\tilde{R}}_i$ and T_i are then concatenated and passed through an embedding network, whose layer dimensions were set to (25, 50, 100) in this work, yielding a local embedding matrix $G_i\in {\mathbb{R}}^{N_i\times M_1}$ , where $N_i$ is the maximum number of neighbors (set to 250 in this work) and M₁ equals 100 which corresponds to the output width of the final embedding layer. The encoded feature $D_i\in {\mathbb{R}}^{M_1\times M_2}$ for atom $i$ is then constructed as:

3

$$D_i={\left(G_i\right)}^T{\tilde{R}}_i{\left({\tilde{R}}_i\right)}^T\stackrel{o}{G_i}$$ where $\stackrel{o}{G_i}$ is a submatrix of $G_i$ to reduce redundancy and T denotes matrix transpose. In our implementation, M₂ was set to 12. D_i is subsequently flattened, concatenated with T_i, and passed into the fitting network for energy fitting. The size of the fitting network was set to (240, 240, 240). Both the embedding and fitting networks include residual connections, as automatically applied by DeePMD-kit when layer widths meet the conditions⁵². All computations were performed using float64 precision.

During training, the loss function L of one single frame is defined by

4

$$L\left(p_e,p_f,p_v\right)=\frac{p_e}{N}\mathrm{\Delta }{E}^2+\frac{p_f}{3N}\sum _i^N\parallel \mathrm{\Delta }{\mathbf{F}}_i{\parallel }^2+\frac{p_v}{9N}\parallel \mathrm{\Delta }\mathbf{V}{\parallel }^2$$ where the N denotes the atom number, Δ denotes the difference between the predicted label and ground truth, E denotes total energy, ΔF_i is the component-wise force difference on atom i, and ΔV is the component-wise difference of the virial tensor. For production MD simulations, the final DeepEMs-25 models were trained for 4,000,000, with the initial learning rate set to 0.001 and decaying every 20,000 steps to finally 3.51 × 10⁻⁸. To assess the robustness of model performance under different data splits, we additionally conducted short-training evaluations. In these trials, models were retrained from scratch using a consistent budget of 1,000,000 steps with learning rate decay every 5000 steps. For each 95/5 random split, training was performed on 95% of the data for 950,000 steps, with decay every 4750 steps, and evaluated on the remaining 5%. These short-trained models were used solely for experimental purposes and not in downstream MD simulations. The similar performance between the random and independent test sets (Supplementary Fig. S13) supports the reliability of the independent test set construction.

Coefficients p_e, p_f, and p_v are dynamic loss weights that gradually evolve during training, starting from (0.1, 1000, 0.2) and converging to (1.0, 1.0, 1.0), respectively. Model compression⁵⁴ was utilized to accelerate the inference procedure.

Descriptor inference and t-SNE

To visualize how the learned descriptors organize atomic environments across energetic materials, a DeepEMs-25 model with PyTorch backend was firstly trained with DeePMD-kit v3⁵⁵ with exactly the same architecture and whole training set as above. Each per-atom descriptor has a fixed dimensionality of 1208, obtained by flattening the $D_i\in {\mathbb{R}}^{100\times 12}$ feature matrix (Eq. (3)) and concatenating an 8-dimensional atomic-type embedding $T_i\in {\mathbb{R}}^8$ . Next, we selected the first 10 frames per material for each of the 20 energetic materials from the test set, computed all per-atom descriptors from these frames, and concatenated them for t-SNE analysis.

All per‑atom latent descriptors for the 20 training materials were first concatenated and paired with their source labels. We then applied principal-component analysis to reduce each high-dimensional descriptor vector to 50 components before computing a two-dimensional t‑SNE with hyperparameters of perplexity set to 80, learning rate set to 200, early exaggeration set to 15 for 250 iterations, and a total 750 iterations.

Efficiency benchmark

We selected DAP-4 for the efficiency benchmark. Supercell expansion was conducted to construct systems with different atom numbers. For DPMD simulations based on Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS), 1000 molecular dynamics steps were executed and the performance metrics were derived by averaging over all steps. In contrast, for the CP2K-based simulations (encompassing both DFT and PM6 methods), a shorter simulation comprising 20 steps was employed to reduce computational expense. To mitigate transient effects associated with the simulation warm-up, only the final 10 steps were averaged to represent the benchmark performance.

Molecular dynamics

MD studies were conducted on the LAMMPS code package^56,57. Thermal decomposition studies were performed in the NVT ensemble with a Nosé-Hoover thermostat at target temperatures. Energy minimization was conducted with specific criteria, followed by 1 ps of NPT pre-equilibrium simulations with an isotropic barostat at 1.0 bar and 300 K. The time step was set to 0.1 fs, with temperature damping parameter set to 0.01 ps for both NPT and NVT simulation, while pressure damping parameter was set to 1000 times the time step (0.1 ps), respectively.

Quantum chemical calculation

(i) The wave functions were calculated using Gaussian 16⁵⁸ with explicitly setting unrestricted PBE0-D3^{49, 50–51,59} functional with def2-TZVP ⁶⁰ basis sets (noted as PBE0-D3/def2-TZVP). The basin analysis for electron density was calculated with Multiwfn 3.8(dev)^41,42 with wavefunctions obtained from Gaussian 16, level PBE0-D3/def2-TZVP. (ii) Zeroth-order symmetry-adapted perturbation based on density fitting (DF-SAPT0) calculations^43,44 were conducted with PSI4⁶¹ with basis set TZVP in considerations of support with heavy atoms like Rb.

Acknowledgements

This study was funded by the National Natural Science Foundation of China (U2341287 and 22488101), Guangzhou Science and Technology Program (2024A04J6499), and Fundamental Research Funds for the Central Universities, Sun Yat-sen University (23lgzy001). The funder played no role in study design, data collection, analysis, interpretation of data, or the writing of this manuscript. The computational resource was supported by the Bohrium Cloud Platform (https://bohrium.dp.tech/) and the National Supercomputing Center in Guangzhou (NSCC-GZ, Tianhe-2). M.-Y.G. thanks Chengqian Zhang and Duo Zhang for fruitful discussions. Language polishing was assisted by DeepSeek-v3 and finalized by human authors.

Author contributions

M.-Y.G. conducted all computational experiments. M.-Y.G. and Y.-F.Y. performed background investigation and data curation. P.C. coordinated high-performance computing resources. W.-X.Z. acquired funding and supervised the entire project.

Data availability

The DeepEMs-25 datasets for training and trained DP models have been uploaded to AIS Square: Datasets–https://aissquare.com/datasets/detail?pageType=datasets&name=DeepEMs-25&id=311. Models – https://aissquare.com/models/detail?pageType=models&name=DeepEMs-25__DPA1-L0&id=312.

Code availability

Codes are openly available in a public repository: https://github.com/SchrodingersCattt/DeepEMs-25-Research. Other major data files are available on request from the authors.

Competing interests

The authors declare no competing interests.

Supplementary information

The online version contains supplementary material available at https://doi.org/10.1038/s41524-025-01739-7.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Hu, Y et al. Chemically driven energetic molecular ferroelectrics. Nat. Commun.; 2021; 12, 1:CAS:528:DC%2BB3MXitFCgtLvF 5696.

2. Chen, C et al. Accurate machine learning models based on small dataset of energetic materials through spatial matrix featurization methods. J. Energy Chem.; 2021; 63, pp. 364-375.1:CAS:528:DC%2BB38XisFOjtbfL

3. Lai, Q et al. Size-matched hydrogen bonded hydroxylammonium frameworks for regulation of energetic materials. Nat. Commun.; 2022; 13, 1:CAS:528:DC%2BB38XivFWksL7L 6937.

4. Li, J et al. Tri-explosophoric groups driven fused energetic heterocycles featuring superior energetic and safety performances outperforms HMX. Nat. Commun.; 2022; 13, 1:CAS:528:DC%2BB38XisFeqtrjO 5697.

5. Zhang, J et al. Synthesis of a high-energy-density material through rapid replacement of crystal water of hydrates. Chem; 2022; 8, pp. 2678-2687.1:CAS:528:DC%2BB38XhvVarsb3M

6. Lee, BH; Larentzos, JP; Brennan, JK; Strachan, A. Graph neural network coarse-grain force field for the molecular crystal RDX. npj Comput. Mater.; 2024; 10, 208.1:CAS:528:DC%2BB2cXitVSrsLvN

7. Wen, L; Wang, Y; Liu, Y. Data-driven combinatorial design of highly energetic materials. Acc. Mater. Res.; 2025; 6, pp. 64-76.1:CAS:528:DC%2BB2cXitlynsL7O

8. Gong, Z et al. High-energy nitrogen rings stabilized by superatomic properties. Adv. Energy Mater.; 2024; 14, 2303446.1:CAS:528:DC%2BB2cXhtFGjsLvJ

9. Chen, S-L et al. Optimizing the oxygen balance by changing the A-site cations in molecular perovskite high-energetic materials. CrystEngComm; 2018; 20, pp. 7458-7463.1:CAS:528:DC%2BC1cXhvFSqtb3O

10. Chen, S-L et al. Molecular perovskite high-energetic materials. Sci. China Mater.; 2018; 61, pp. 1123-1128.1:CAS:528:DC%2BC1MXhvV2msb0%3D

11. Shang, Y et al. Metal-free molecular perovskite high-energetic materials. Cryst. Growth Des.; 2020; 20, pp. 1891-1897.1:CAS:528:DC%2BB3cXivVaktrk%3D

12. Zhang, W-X; Chen, S-L; Shang, Y; Yu, Z-H; Chen, X-M. Molecular perovskites as a new platform for designing advanced multi-component energetic crystals. Energ. Mater. Front.; 2020; 1, pp. 123-135.

13. Feng, X-J; Zhang, K; Shang, Y; Pan, W. Thermal decomposition study of perchlorate-based metal-free molecular perovskite DAP-4 mixed with ammonium perchlorate. Case Stud. Therm. Eng.; 2022; 34, 102013.

14. Wang, X; Wang, Q; Wang, X; Wang, J; Zhang, G. Thermal pyrolysis properties and security performance of molecular perovskite energetic crystal (C₆N₂H₁₄)(NH₄)(ClO₄)₃ (DAP-4). J. Anal. Appl. Pyrolysis; 2023; 176, 106268.1:CAS:528:DC%2BB3sXis1alurvO

15. Wang, X et al. Thermal properties and decomposition of perovskite energetic materials (C₆H₁₄N₂)NH₄(ClO₄)₃. Thermochim. Acta; 2024; 742, 1:CAS:528:DC%2BB2cXit1Ghs7%2FJ 179878.

16. Zhou, J et al. Comparative thermal research on energetic molecular perovskite structures. Molecules; 2022; 27, 805.1:CAS:528:DC%2BB38XktFWrtb0%3D

17. Deng, P; Wang, H; Yang, X; Ren, H; Jiao, Q. Thermal decomposition and combustion performance of high-energy ammonium perchlorate-based molecular perovskite. J. Alloy. Compd.; 2020; 827, 154257.1:CAS:528:DC%2BB3cXktVGmurw%3D

18. Zhai, P et al. Thermal decomposition of ammonium perchlorate-based molecular perovskite from TG-DSC-FTIR-MS and ab initio molecular dynamics. RSC Adv.; 2021; 11, pp. 16388-16395.1:CAS:528:DC%2BB3MXhtVWqu7jE

19. Senftle, TP et al. The ReaxFF reactive force-field: development, applications and future directions. npj Comput. Mater.; 2016; 2, 15011.1:CAS:528:DC%2BC2sXlslantL4%3D

20. Van Duin, AC; Dasgupta, S; Lorant, F; Goddard, WA. ReaxFF: a reactive force field for hydrocarbons. J. Phys. Chem. A; 2001; 105, pp. 9396-9409.

21. Schreiner, M; Bhowmik, A; Vegge, T; Busk, J; Winther, O. Transition1x-a dataset for building generalizable reactive machine learning potentials. Sci. Data; 2022; 9, 779.

22. Zhang, S; Zubatyuk, R; Yang, Y; Roitberg, A; Isayev, O. ANI-1xBB: an ANI-based reactive potential for small organic molecules. J. Chem. Theory Comput.; 2025; 21, pp. 4365-4374.1:CAS:528:DC%2BB2MXpsF2jtbg%3D

23. Levine, D. S. et al. The open molecules 2025 (OMol25) dataset, evaluations, and models. Preprint at https://arxiv.org/abs/2505.08762 (2025).

24. Wen, M. et al. A general neural network potential for energetic materials with C, H, N, and O elements. Preprint at https://arxiv.org/abs/2503.01932 (2025).

25. Wang, H; Zhang, L; Han, J; Weinan, E. DeePMD-kit: a deep learning package for many-body potential energy representation and molecular dynamics. Comput. Phys. Commun.; 2018; 228, pp. 178-184.1:CAS:528:DC%2BC1cXmt1ajsrk%3D

26. Zhang, L; Han, J; Wang, H; Car, REW. Deep potential molecular dynamics: a scalable model with the accuracy of quantum mechanics. Phys. Rev. Lett.; 2018; 120, 143001.1:CAS:528:DC%2BC1MXltFSksrg%3D

27. Zhang, L., Han, J., Wang, H., Saidi, W. & Car, R. End-to-end symmetry preserving inter-atomic potential energy model for finite and extended systems. In Advances in Neural Information Processing Systems, 31 (NeurIPS, 2018).

28. Sumpter, BG; Noid, DW. Potential energy surfaces for macromolecules. a neural network technique. Chem. Phys. Lett.; 1992; 192, pp. 455-462.1:CAS:528:DyaK38XktVyru7s%3D

29. Behler, J; Parrinello, M. Generalized neural-network representation of high-dimensional potential-energy surfaces. Phys. Rev. Lett.; 2007; 98, 146401.

30. Zhang, D et al. Pretraining of attention-based deep learning potential model for molecular simulation. npj Comput. Mater.; 2024; 10, 94.

31. Smith, JS; Nebgen, B; Lubbers, N; Isayev, O; Roitberg, AE. Less is more: sampling chemical space with active learning. J. Chem. Phys.; 2018; 148, 241733.

32. Zhang, L; Lin, D-Y; Wang, H; Car, R; E, W. Active learning of uniformly accurate interatomic potentials for materials simulation. Phys. Rev. Mater.; 2019; 3, 023804.1:CAS:528:DC%2BC1MXhtFGrsrnF

33. Zhang, Y et al. DP-GEN: a concurrent learning platform for the generation of reliable deep learning based potential energy models. Comput. Phys. Commun.; 2020; 253, 107206.1:CAS:528:DC%2BB3cXjvFWjsrk%3D

34. O’Boyle, NM et al. Open Babel: an open chemical toolbox. J. Cheminformatics; 2011; 3, 33.

35. Zeng, J et al. ReacNetGenerator: an automatic reaction network generator for reactive molecular dynamics simulations. Phys. Chem. Chem. Phys.; 2020; 22, pp. 683-691.1:CAS:528:DC%2BC1MXit1artLnM

36. Wu, Y; Sun, H; Wu, L; Deetz, JD. Extracting the mechanisms and kinetic models of complex reactions from atomistic simulation data. J. Comput. Chem.; 2019; 40, pp. 1586-1592.1:CAS:528:DC%2BC1MXjs1Gmurw%3D

37. Du, P; Zhu, H; Braun, A; Yelon, A; Chen, Q. Entropy and isokinetic temperature in fast ion transport. Adv. Sci.; 2024; 11, 2305065.1:CAS:528:DC%2BB3sXit1yrtbvL

38. Takamure, N et al. Thermodynamic interpretation of the Meyer-Neldel rule explains temperature dependence of ion diffusion in silicate glass. Phys. Rev. Lett.; 2022; 129, 1:CAS:528:DC%2BB38XivVWlsb7L 175901.

39. Shashidhar, MS; Krishnaswamy, S. Intermolecular acyl-transfer reactions in molecular crystals. Acc. Chem. Res.; 2019; 52, pp. 437-446.1:CAS:528:DC%2BC1MXitVSitro%3D

40. Bader, RF. Atoms in molecules. Acc. Chem. Res.; 1985; 18, pp. 9-15.1:CAS:528:DyaL2MXmtFGgsA%3D%3D

41. Lu, T; Chen, F. Multiwfn: a multifunctional wavefunction analyzer. J. Comput. Chem.; 2012; 33, pp. 580-592.

42. Lu, T. A comprehensive electron wavefunction analysis toolbox for chemists, Multiwfn. J. Chem. Phys.; 2024; 161, 1:CAS:528:DC%2BB2cXhvVyiur%2FE 082503.

43. Hohenstein, EG; Sherrill, CD. Density fitting and Cholesky decomposition approximations in symmetry-adapted perturbation theory: Implementation and application to probe the nature of π-π interactions in linear acenes. J. Chem. Phys.; 2010; 132, 184111.

44. Hohenstein, EG; Parrish, RM; Sherrill, CD; Turney, JM; Schaefer, HF. Large-scale symmetry-adapted perturbation theory computations via density fitting and Laplace transformation techniques: investigating the fundamental forces of DNA-intercalator interactions. J. Chem. Phys.; 2011; 135, 174107.

45. Anderson, PW. More Is Different: Broken symmetry and the nature of the hierarchical structure of science. Science; 1972; 177, pp. 393-396.1:CAS:528:DyaE38XltVGlu7s%3D

46. Kühne, TD et al. CP2K: an electronic structure and molecular dynamics software package-Quickstep: Efficient and accurate electronic structure calculations. J. Chem. Phys.; 2020; 152, 194103.

47. VandeVondele, J; Hutter, J. An efficient orbital transformation method for electronic structure calculations. J. Chem. Phys.; 2003; 118, pp. 4365-4369.1:CAS:528:DC%2BD3sXhtlygs7g%3D

48. Weber, V; VandeVondele, J; Hutter, J; Niklasson, A. Direct energy functional minimization under orthogonality constraints. J. Chem. Phys.; 2008; 128, 084113.

49. Perdew, JP; Burke, K; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett.; 1996; 77, 3865.1:CAS:528:DyaK28XmsVCgsbs%3D

50. Grimme, S; Antony, J; Ehrlich, S; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H–Pu. J. Chem. Phys.; 2010; 132, 154104.

51. Grimme, S; Ehrlich, S; Goerigk, L. Effect of the damping function in dispersion corrected density functional theory. J. Comput. Chem.; 2011; 32, pp. 1456-1465.1:CAS:528:DC%2BC3MXjsF2isL0%3D

52. Zeng, J et al. DeePMD-kit v2: A software package for deep potential models. J. Chem. Phys.; 2023; 159, 054801.1:CAS:528:DC%2BB3sXhs1Sksb7L

53. Lu, D et al. 86 PFLOPS deep potential molecular dynamics simulation of 100 million atoms with ab initio accuracy. Case Stud. Therm. Eng.; 2021; 259, 1:CAS:528:DC%2BB3cXitVKgt7vI 107624.

54. Lu, D et al. DP compress: A model compression scheme for generating efficient deep potential models. J. Chem. Theory Comput.; 2022; 18, pp. 5559-5567.1:CAS:528:DC%2BB38XitVCgt7rF

55. Zeng, J et al. DeePMD-kit v3: a multiple-backend framework for machine learning potentials. J. Chem. Theory Comput.; 2025; 21, pp. 4375-4385.1:CAS:528:DC%2BB2MXhtVChu7vM

56. Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys.; 1995; 117, pp. 1-19.1:CAS:528:DyaK2MXlt1ejs7Y%3D

57. Thompson, AP et al. LAMMPS-a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. Comput. Phys. Commun.; 2022; 271, 1:CAS:528:DC%2BB3MXitlSrsb7O 108171.

58. Frisch, M. et al. Gaussian 16 Revision C. 01. 2016, 421 (Gaussian Inc. Wallingford CT, 2016).

59. Adamo, C; Barone, V. Toward reliable density functional methods without adjustable parameters: the PBE0 model. J. Chem. Phys.; 1999; 110, pp. 6158-6170.1:CAS:528:DyaK1MXitVCmt7Y%3D

60. Weigend, F; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. Phys. Chem. Chem. Phys.; 2005; 7, pp. 3297-3305.1:CAS:528:DC%2BD2MXpsFWgu7o%3D

61. Turney, JM et al. Psi4: an open-source ab initio electronic structure program. WIREs Comput. Mol. Sci.; 2012; 2, pp. 556-565.1:CAS:528:DC%2BC38XhtlWiu7bF

62. Shannon, RD. Revised effective ionic radii and systematic studies of interatomic distances in halides and chalcogenides. Found. Crystallogr.; 1976; 32, pp. 751-767.

63. Sidey, V. On the effective ionic radii for ammonium. Struct. Sci.; 2016; 72, pp. 626-633.1:CAS:528:DC%2BC28Xht1ymtr7E