All-nitrogen cages and molecular crystals: Topological rules, stability, and pyrolysis paths

We have combined ab initio molecular dynamics with the intrinsic reaction coordinate in order to investigate the mechanisms of stability and pyrolysis of N$_{4}$-- N$_{120}$ fullerene-like nitrogen cages. The stability of the cages was evaluated in terms of the activation barriers and the activation Gibbs energies of their thermal-induced breaking. We found that binding energies, bond lengths, and quantum-mechanical descriptors failed to predict the stability of the cages. However, we derived a simple topological rule that adjacent hexagons on the cage surface resulted in its instability. For this reason, the number of stable nitrogen cages is significantly restricted in comparison with their carbon counterparts. As a rule, smaller clusters are more stable, whereas earlier proposed rather large cages collapse at room temperature. The most stable all-nitrogen cages are N$_{4}$ and N$_{6}$ clusters, which can form the van-der-Waals crystals with the densities of 1.23 and 1.36 g/cm$^{3}$, respectively. Examination of their band structures and densities of electronic states shows that they are both insulators. Their power and sensitivity are not inferior to the modern advanced high-energy nanosystems.


Introduction
The formation of the N2 molecule with the triple N≡N bond from two isolated nitrogen atoms results in the release of a large amount of energy (~ 10 eV [1]). Thus, nitrogen is considered as the primary component of the most high-energy compounds. Such compounds usually contain also carbon atoms that provide stability of the whole molecule framework decorated by oxygencontaining groups or other oxidants. RDX, HMX, and CL-20 are well-known examples of traditional high-energy structures with a carbon-nitrogen frame. In recent years, more advanced energy materials have been proposed with a similar carbon-nitrogen architecture [2; 3; 4].
The recent synthesis of non-molecular nitrogen [5; 6] has renewed research interest in allnitrogen structures without carbon additives. Compared to carbon-nitrogen compounds, pristine nitrogen systems are more powerful and safer for the environment. In particular, closed fullnitrogen Nn cages with three-coordinated atoms and single N-N bonds are very attractive and environmentally friendly oxygen-free fuels with an extremely high energy capacity. The transformation of single N-N bonds into the triple bond is a very exothermic process. The breakdown of Nn cages into N2 molecules gives more than 50 kcal/mol per nitrogen atom [7].
However, in practice, Nn cages possess low stability and have not yet been synthesized and measured. Unlike carbon fullerenes, fully nitrogen cages do not retain their structure under ambient conditions, since they are energetically unfeasible compared to N2 molecules. In addition, the activation energies of their decay processes are rather low. Thus, most Nn cages can be prepared and stabilized only at extremely high pressures, similar to recently obtained forms of nonmolecular nitrogen [5; 6]. The introduction of carbon atoms, oxygen, or other atoms into the nitrogen cages leads to more stable structures [3; 4; 8], which, however, are not as efficient and environmentally friendly as all-nitrogen hypothetical analogs [7; 9].
Low kinetic stability can be considered as the Achilles heel of all previously studied systems, even if they were true energy minima in the configuration space.
In the presented paper, we carry out a systematic study of the structure and stability of Nn fullerene-like cages. Our goal is to identify the general relationships between topology and stability of nitrogen buckyballs. Unlike previous calculations, we focused on the pyrolysis mechanisms and the corresponding activation barriers, rather than on the relative formation energies of the structures under consideration. We identified the most stable Nn cages, explained the reasons for their stability, and investigated the possibilities of their crystallization.

Computational Details
Finite-size nitrogen cages were simulated in the frame of the density functional theory with the B3LYP functional [23; 24] and 6-311G(d,p) basis set [25]. GAMESS software [26]   where  and  are real eigenfrequencies at stable and transition states, respectively. Note that the imaginary eigenfrequency of the transition state is not included in the multiplication.
Degeneracy factor g is equal to the number of equivalent decay paths (g > 1 for high-symmetry Nn cages).
Ab initio molecular dynamics was carried out using a high-performance GPU-based TeraChem package [28; 29; 30; 31] with the same B3LYP/6-311G(d,p) level of theory. The initial atomic displacements and velocities were generated in accordance with the desired temperature.
During the simulation, the given temperature was maintained by the Langevin thermostat [32].
The time step for all simulations was 0.1 fs, which was sufficient to take into account all molecular vibrations correctly. Systems' dynamical trajectories were visualized using the ChemCraft program [33].
To study geometry and electronic structure properties of crystals that consist of the nonmolecular nitrogen systems, we used another density functional theory approach and its implementation in the QUANTUM Espresso 6.5 program package [34; 35]. The generalized gradient approximation (GGA) in the Perdew-Burke-Ernzerhof (PBE) functional form for the exchange-correlation energy [36], and the projector-augmented-wave (PAW) method for the electron-ion interaction [37; 38] was used to perform the calculations. The value of cut-off energy for the plane-wave basis set was chosen as 120 Ry (1632 eV). The weak van der Waals interactions between the non-covalently bound nitrogen atoms are taken into account by using the D3 Grimme dispersion corrections [39]. The geometry optimization of the nitrogen crystal was performed without any symmetry constraints until the Hellman-Feynman forces acting on the atoms became smaller than 10 -6 hartree/bohr. Supercell parameters were also optimized. The Brillouin zone integrations were performed by using the Monkhorst-Pack k-point sampling scheme [40] with the 8×8×8 mesh grid. For the non-self-consistent field calculations, the k-point grid size of 24×24×24 was used. The Methfessel-Paxton smearing [41] was used for the geometry relaxation with the smearing width of 0.02 eV. Still, for the calculations of the electronic density of states, the Böchl tetrahedron method [42] was employed. The electronic structure properties were elucidated through the analysis of the sample band structure and its electronic density of states.

A. The first insight on the bicyclic hydro-nitrogen molecules
Hypothetical bicyclic hydro-nitrogen molecules are the smallest possible systems containing adjacent nitrogen rings with the single N-N bonds. We constructed N10H8, N9H7, N8H6, N7H5, and N6H4 naphthalene-like molecules with hexagon/hexagon, hexagon/pentagon, pentagon/pentagon, pentagon/square, and square/square interfaces, respectively (see Fig. 1a).
Unlike their carbon counterparts, these molecules did not remain flat during geometry optimization.
We found that the hexagon/pentagon and the pentagon/pentagon were the only stable configurations with positive eigenfrequencies. Other structures are not true energy minima.
Potential energy paths corresponding to the breaking of joint hexagon/pentagon and pentagon/pentagon bonds are shown in Fig. 1b and 1c (the energy barriers are 0.58 and 0.47 eV, respectively). Therefore, the pentagon/pentagon interface seems to be the most suitable. The hexagon/pentagon interface possesses lower stability, while the adjacent hexagons lead to the complete instability. Note that the carbon structures exhibit the opposite behavior, avoiding adjacent pentagons according to Kroto's isolated pentagon rule [43].
However, bicyclic molecules can only give qualitative insight into the stability of allnitrogen cages. They do not reproduce local curvature effects, which is especially crucial for small cages, and other effects such as three-dimensional aromaticity. Moreover, hydrogen passivation distorts the charge distribution in the nitrogen rings. Thus, the conclusions made from the analysis of bicyclic nitrogen-containing molecules should be verified on more suitable all-nitrogen systems.

B. Nitrogen cages with the adjacent hexagons
As a next step, we constructed nitrogen analogs of low-energy isomers of carbon fullerenes from C20 to C120, created using the FULLERENE 4.5 code [44]. All possible isomers C20 -C32, as well as all low-energy isomers C34 -C120, included in the Tomanek database [45] were used as the initial sample geometry for all-nitrogen cages. Note that the typical length of a single N-N bond  Fig. 2a. The extremely small diameters of these tubes provide stability due to the curvature effects and concave hexagon shapes. Note that the large diameter Nn tubular structures are unstable [46]. We test the stability of the tubular systems shown in Fig.2 using ab initio molecular dynamics (AIMD) at T = 700 K for one picosecond. We observed their decomposition into N2 molecules, and short Nn chains (n = 3 ÷ 7) started with breaking the adjacent sides of the neighboring hexagons. We conclude that their stability is too low even for their identification. Such nitrogen nanostructures can hardly be used for any practical applications under ambient conditions. Our calculations confirm that, unlike carbon, nitrogen avoids the formation of adjacent hexagonal rings. Note that at the moment, compounds containing adjacent N6 rings have not yet been experimentally synthesized.

C. Nitrogen cages with the isolated hexagons
A large number of hexagons on the fullerene surface leads to a small curvature and large size of the system. To avoid adjacent pentagons, carbon clusters Cn should be sufficiently large (n ≥ 60). In contrast, the hexagons on the surface of nitrogen clusters should be separated by pentagons or smaller polygons, and therefore the fracture of the hexagons should be low. Thus, stable nitrogen clusters Nn should have a significant curvature and small size (n ≤ 24). Inspired by small carbon fullerenes and fullerene-like cages (including tetrahedranes, prismanes, non-classical fullerenes, etc.), we have constructed the corresponding nitrogen structures. All the nitrogen cages satisfying the isolated hexagon rule that are true local energy minima are presented in Fig. 3. The kinetic stability of these systems was investigated using AIMD at T = 700 K for one picosecond.
Only the five smallest nitrogen clusters, labeled as (a) -(e) in Fig. 3, conserved their identity during the AIMD simulation (their pyrolysis mechanisms will be described in more detail in the next subsection). Larger clusters are decomposed. Therefore, they should be considered as low-stable systems. Note that N-N bonds in all truly stable and low-stable cages are longer than 1. To get a complete understanding of the origin of the stability of nitrogen cages, we analyzed the mechanisms of their pyrolysis. We adopted the reaction mechanism from the ab initio Thermodynamic descriptors (activation enthalpies and Gibbs activation energies) at various temperatures are collected in Table 1. These descriptors weakly depend on temperature for all clusters, since the initial and transition configurations are not so different, and their vibrational energies and entropies slightly differ from each other.  [47] (for example, inside porous materials or carbon cages [48]).
Due to the compact shape of small nitrogen clusters, spatial confinement can increase their decomposition activation barrier by 0.5 eV or more [48] that is sufficient to stabilize them at room temperature but hardly useful at higher temperatures.

E. Van der Waals nitrogen crystals
Recently, some crystal structures have appeared based on the chains N6 [49] and N8 [50], as well as other single-bonded all-nitrogen systems [51; 52]. According to the above calculations, the N4 and N6 clusters are very stable and can be considered as an alternative to nitrogen chains.
Their activation barriers are even higher than the corresponding values for experimentally observed highly strained hydrocarbons (for example, tetrahedrane C4H4 [53] or polymethylcubanes C8Hq(CH3)8-q [54]). Thus, N4 and N6 are the most suitable candidates for building blocks for the high-energy-density crystals.
There are many algorithms for predicting the crystal structure. The most famous of them are USPEX [55] and CALYPSO [56]. However, they seem rather excessive for the construction of relatively simple molecular crystals in which the N4 and N6 cages conserve their shape. Isolated nitrogen cages bind in such crystals due to the weak Van der Waals interaction. Instead, we constructed three initial lattices for each crystal (sc, bcc, and fcc) and optimized atomic positions along with lattice vectors without any symmetry constraints. Then the crystals with the lowest potential energy were considered. Their crystallographic parameters are summarized in Table 2.
Their band structures and electron densities of states are shown in Fig. 6. Coordinates of the atomic positions are available in Supplementary Materials.
Assuming the detonation of these crystals with the release of gaseous N2, we estimated the heat of detonation Q as the difference in energies between crystals and nitrogen molecules. In addition, the characteristics of detonation (detonation velocity D, adiabatic exponent Γ, detonation pressure P) can be estimated from Q and crystal density p in accordance with the Xiong empirical rules [57] that, in our case, have the following forms:  Table 2), which is two times higher than the corresponding value for ε-CL-20 (12 cm [61]). Note that this is only a qualitative estimation of impact sensitivity.
However, this confirms the potential applicability of the nitrogen crystals under consideration.

CONCLUSION
In this paper, we analyzed the structure and stability of all-nitrogen fullerene-like clusters and molecular crystals based on them. According to the formulated topological rule, confirmed by density functional theory calculations, only small nitrogen cages with single N-N bonds are stable.
In terms of stability, they can compete with actively studied nitrogen chains. An increase in cage size always leads to a sharp decrease in its stability. Thus, the range of nitrogen cages is significantly limited. Larger structures can be formed from small building blocks due to noncovalent interactions.
However, our conclusions are valid only for convex cages. We admit that the concave shape of the nitrogen structure can contribute to its stability. However, the curvature of such a hypothetical construction must change its shape many times in order to avoid large convex regions.
The study of such astralene-like structures is beyond the scope of this study. We believe that the topological rules presented here will facilitate the further search for new clusters and crystals of nitrogen.
Like strained hydrocarbon frameworks, small nitrogen clusters are stable enough to form molecular crystals. The performed analysis of the two most promising all-nitrogen structures confirmed their applicability as competitive, high-energy materials.

ACKNOWLEDGMENTS
The reported study was funded by RFBR according to the research project No. 18-32-