Molecular nitrogen (N₂) is the most abundant component of Earth’s atmosphere and one of the most stable molecules owing to its strong N≡N triple bond. There is a large energy difference between single/double nitrogen bonds (160 and 418 kJ/mol for N–N and N=N, respectively) and the N≡N triple bond (954 kJ/mol).1 Therefore, polynitrogen composed of mixed single and double bonds can store a substantial amount of chemical energy. The stored energy will be released when decomposition of the polymeric nitrogen is triggered, producing pollution-free N₂ molecules, in an oxygen-free process. However, the high bond energy of the N₂ molecule makes it highly unreactive. It is difficult to synthesize polymeric nitrogen from N₂ molecules. To break nitrogen triple bonds and obtain nitrogen polymers, extreme conditions such as high pressures and high temperatures become necessary to overcome the high energy barriers.2,3 For example, pressure can be used to break the nitrogen bonds, with the formation of some unexpected nitrogen compounds.4 In some nitrogen-rich cases, pressure can enhance the effect of electronic delocalization, thus helping to dissociate molecules and form polymeric structures.

Tremendous efforts have been made to obtain novel polymeric nitrogen structures under high pressures in both theoretical and experimental studies.2,3,5–9 More than a decade after its theoretical prediction by Mailhiot et al.,5 an experimental breakthrough was eventually made by Eremets et al. in 2004 in their synthesis of a cubic-gauche nitrogen structure (cg-N).2 This cg-N contains pure N–N single bonds and thus possesses an extremely high energy density, which is reported to be five times higher than that of the current most powerful energetic materials.5 Since then, several polymeric nitrogen structures have been synthesized: hexagonal layered polymeric nitrogen (HLP-N), layered polymeric nitrogen (LP-N), and a black phosphorus structure (BP-N).10–13 However, the pressure–temperature conditions required for these polymeric nitrogen syntheses are very high, exceeding 110 GPa and 2000 K. Such extreme conditions cause great practical difficulties, and therefore the search for potential high-energy-density materials (HEDMs) that can be synthesized under moderate pressures is still an important task.

Recently, counter-ions and nitrogen fragments have been used to construct nitrogen-bearing compounds involving ionic bonding, including N₅ and N₆ rings, nitrogen azide, and polymeric N₄ chains (poly-N42−).7,14–28 In these cases, the ionic bonds can enhance the kinetic stability of the polynitride components and reduce the pressure required for synthesis. Among these nitrogen fragments, poly-N42− has become of particular importance. Previous theoretical works have predicted the existence of several poly-N42− compounds, including BeN₄, MgN₄, CaN₄, CdN₄, and FeN₄.29–33 Of these, MgN₄, FeN₄, TaN₄, and BeN₄ have recently been synthesized.34–37 Typical alkaline-earth nitrides such as the BeN₄, MgN₄, and CaN₄ compounds,29,30,32 have similar N–N bond lengths in the range 1.32–1.35 Å. However, the poly-N42− chains in FeN₄34 possess significantly different N–N bond lengths, between 1.29 and 1.43 Å. Nevertheless, these poly-N42− chains are constructed with N–N single bonds and N=N double bonds, matching the features of HEDM nitrides. Additionally, the N2− species is found to coexist with the poly-N42− chain in some high-energy-density metal nitrides synthesized by Bykov et al.38,39 The cations in these structures are mainly divalent metals. We envision that the nitrogen content of metal nitrides can be increased by the use of metal ions in higher valence states, thereby achieving higher energy densities.

In this work, we focus on the metals Al, Ga, Y and Ti, which can provide three or four valence electrons. We predict that the poly-N42− chain can occur in four new metal nitride crystals MN_x with different symmetries: P2₁-AlN₆, P2₁-GaN₆, P-1-YN₆, and P4/mnc-TiN₈. These are predicted to be energetically stable at moderate pressures. Molecular dynamics (MD) simulations and phonon spectrum calculations suggest that these MN_x are mechanically stable at nonzero temperatures and high pressures. More interestingly, we find coexistence of the aforementioned two types of poly-N42− chains, appearing in our predicted P2₁ phase of the AlN₆ and GaN₆ compounds. Besides, we estimate that these MN_x have excellent detonation properties and can served as potential HEDMs.

We carry out a crystal structure search using MAGUS (“machine learning and graph theory assisted universal structure searcher”),40,41 which is accelerated through the use of graph theory42 and machine learning potentials. This method has been successfully applied to many systems.40,43–48 To generate new structures containing the desired N42− units, molecular-based unit searches are conducted at 50 GPa. We constrain four MN_x formula units into one cell to explore more nitrogen-containing configurations. The first-principles calculations are performed using the Vienna Ab initio Simulation Package (VASP).49 Density-functional theory (DFT) is implemented within the projector-augmented-wave (PAW) approach.50 We use the Perdew–Burke–Erzernhof functional based on the generalized gradient approximation51,52 to describe the exchange–correlation interaction. Van der Waals (vdW) effects are taken into account using the DFT-D3 method in VASP.53,54 The valence electrons occupying the shell states are 3s²3p¹, 4s²4p¹, 4s²4p⁶5s¹4d², and 3p⁶4s²3d² for Al, Ga, Y, and Ti, respectively. The energy cutoff of the plane-wave basis is set at 1050 eV. The Brillouin zone is separated with k-meshes of 2π × 0.03 Å⁻¹ to guarantee convergence of energy. We use the PHONOPY code55 to calculate the phonon spectrum for a 2 × 2 × 2 supercell. Projected crystal orbital Hamilton population (pCOHP) calculations are conducted using the LOBSTER code.56,57 Noncovalent interactions are also taken into consideration by calculating the reduced density gradient (RDG) using the Critic2 code.58,59 The structures are visualized by VESTA software.60 We analyze the electronic properties in the context of molecular orbital (MO) theory.61

Using the crystal structure search method described above, we predict several metal nitrides structures, as shown in Figs. 1(a)–1(c). Diverse channel frames are constructed from different coordinate bonds between metal atoms and poly-N42− chains [Figs. 1(d)–1(f)]. These channel-like structures may capture small molecules, such as N₂ or H₂O, all of which stabilize these metal nitrides.17,18,39 For AlN₆ and GaN₆ in the P2₁ phase, one metal atom is found to be coordinated with three or four chains of poly-N42− [Fig. 1(d)]. As shown in Fig. 1(e), one metal atom in the P-1-YN₆ structure is sevenfold-coordinated with nitrogen atoms. Because N₂ double balls coexist with poly-N42− chains in the P-1-YN₆ structure, we can describe the YN₆ within a unit cell as Y4(N2)2(N4)5. For P4/mnc-TiN₈, the poly-N42− chain is fourfold-coordinated around one Ti atom [Fig. 1(f)].

We calculate the relative enthalpy of MN_x (M = Al, Ga, Y, and Ti) as shown in Fig. 2. By taking the enthalpy difference of MN_x relative to the ground-state metal nitrides (MN) plus bulk nitrogen phase at these pressures, ΔH = MN_x − MN − (x − 1)N, we can confirm the energetic stabilities of MN_x compounds under pressure. Here, the Pa-3 (α) and P4₁2₁2 phases of the N₂ molecular crystal structures are selected for the relative-enthalpy calculations at 0 GPa and at higher pressures, respectively.65Fm-3m-MN is chosen as the ground state for the Al–N, Y–N and Ti–N systems. The ground state for the GaN structure in the Ga–N system favors the P6₃/mmc phase, which tends to transform into the Fm-3m phase above 32.5 GPa.62 We find that AlN₆, GaN₆, and TiN₈ are enthalpically stable above pressures of around 40.8, 39.3, and 40.1 GPa, respectively. However, YN₆ turns out to be stable above a much lower external pressure (around 21.7 GPa). Therefore, we suggest that MN_x might be synthesized at moderate pressures around 40 GPa or lower. In addition, no imaginary frequency of the phonon spectrum has been found in the Brillouin zones of these MN_x structures, indicating their dynamical stability, as shown in Fig. S1 in the supplementary material. Among them, in particular, P4/mnc-TiN₈ is quenchable to 0 GPa.

To confirm the thermal stabilities of the newly predicted MN_x, we conduct ab initio molecular dynamics (AIMD) simulations for the MN_x (2 × 2 × 2) supercells at high temperatures for 12 ps with an interval of 1 fs. The supercells contain 224 atoms for the AlN₆, GaN₆, and YN₆ structures and 144 atoms for TiN₈. The AIMD calculations are performed within the canonical ensemble for AlN₆ and GaN₆ at 20 GPa, for YN₆ at 40 GPa, and for TiN₈ at 30 GPa. The resulting radial distribution functions (RDFs) are presented in Fig. 3 and in Fig. S2 in the supplementary material. The sharp peak around 1.30 Å in Fig. 3 represents the N–N bonds within N₄ units in the balanced structures. It matches well with the N–N bond length in the original structures (the vertical dashed line). Therefore, the N₄ units remain intact to at least 1200 K for AlN₆, GaN₆, and TiN₈, and 900 K for YN₆. The insets show the stable averaged structures at the corresponding temperatures. In addition, the RDFs of metal–nitrogen (M–N) pairs and metal–metal (M–M) pairs (Fig. S2 in the supplementary material) suggest that the M–N and M–M pairs remain stable up to the corresponding temperatures. At higher temperatures, the MN_x frameworks begin to decompose. We use the mean-squared displacement (MSD), which represents the overall particle deviation, to measure whether the MN_x break down. The MSD function is defined by MSD(t) = u〈|r_i(t) − r_i(0)|²〉, where r_i(t) and r_i(0) are the positions of atom i at time t and the initial time, respectively. Figure S3 in the supplementary material presents the MSDs for AlN₆ at 3000 K, GaN₆ at 2500 K, YN₆ at 1200 K, and TiN₈ at 2500 K. The increases in MSDs mean that the MN_x exhibit diffusive tendencies. Also, as shown by Fig. S4 in the supplementary material, the structures of MN_x at 12 ps exhibit large deviations from the original structures. Some N–N bonds in poly-N42− chains are stretched or broken. Depending on the degree of decomposition, the poly-N42− chains become separate N₄, N₃, or N₂ units, and so we can deduce that MN_x will break down at the corresponding temperatures.

The calculated density of states (DOS) exhibits metallic features for all MN_x structures, as shown in Figs. 4(a)–4(d). The main contribution to the DOS at the Fermi energy is from the N-2p orbital for the AlN₆ and GaN₆ structures, while the d orbitals of the transition metals Y and Ti make additional contributions to the DOS around the Fermi energy for YN₆, and TiN₈. We carry out a projected crystal orbital Hamilton population (pCOHP) calculation to distinguish the bond features in the predicted MN_x nitrides. Using minus pCOHP (−pCOHP), we can partition the bond energy into bonding states with positive values and antibonding states with negative values. Figures 4(e)–4(h) reveal similar bond features for these four structures. For N–N bonds, although most of the states below the Fermi level can be attributed to bonding states, there are still some antibonding states existing just below the Fermi level. For M–N bonds, however, almost all the states below the Fermi level are bonding states. The minus integral pCOHP (−IpCOHP)66 values at the Fermi level, which represent the pairwise interatomic interaction strength, are listed in Table I. The −IpCOHP values of N–N bonds are at least twice those of M–N bonds, which suggests that the N–N bonds have greater pairwise interatomic interaction strengths than the M–N bonds.

To explain the electronic features of MN_x, we apply MO theory61 to analyze the electron orbitals and describe the coordination of the metal atoms. For the isostructural AlN₆ and GaN₆, two types of poly-N42− chains (types A and B) are sketched in Figs. 5(a) and 5(b), respectively. Each N atom in type A [Fig. 5(a)] is surrounded by three atoms in a nearly planar triangular geometry, suggesting sp² hybridization. Of a total of 22 electrons in one N42− unit, 16 are found to occupy eight sp² hybrid orbitals. The extra four p orbitals of four nitrogen atoms form four π orbitals. The remaining electrons occupying π orbitals form delocalized π bonds, which is the origin of metallicity. The equally distributed π band also explains why the N–N bond lengths (1.32–1.33 Å) are between those of an N–N single bond (1.45 Å) and an N=N double bond (1.20 Å).29 As shown in Fig. 5(b), for the type B poly-N42− chain, the two N atoms colored red possess tetrahedral coordination including lone pairs, indicating their sp³ hybridization, while the other two N atoms are sp²-hybridized, similar to those in Fig. 5(a), and therefore the extra two p atomic orbitals of these two N atoms form one bonding orbital π and one antibonding orbital π*. The remaining two electrons occupying the bonding π orbital form a π bond, distributed over the sp²-hybridized N atoms. This can explain why the length of the bond between the two sp³-hybridized N atoms is 1.40 Å, which is close to that of an N–N single bond (1.45 Å), but the length of the bond connecting the sp²-hybridized N atoms is 1.30 Å, which is between those of single and double bonds. Therefore, the N–N bonds in the type A chain are all hybridized bonds, while those in the type B chain are half-single, half-hybridized. In terms of bond energy, the type B chain may possess a higher energy capability in principle. However, YN₆ and TiN₈ contain only type A poly-N42− chains.

When the electron density ρ is too low to allow analysis of the results for the electron localization function, we can identify the noncovalent interactions (NCIs) by calculating the reduced density gradient (RDG) s, which is defined ass=12(3π2)1/3|∇ρ|ρ4/3.As shown in Fig. 6, a plot of s vs ρ multiplied by the sign of the second Hessian eigenvalue of ∇²ρ, i.e., sign(λ₂),59,67 reveals the basic partition of the intramolecular interactions, including weak NCIs. Here, sign(λ₂) identifies the attractive weak interactions (negative, far left), repulsive weak interactions (positive, far right), and van der Waals interactions (around zero). The obvious spikes at low ρ, with much faster decrease of the gradient of ρ, indicate the existence of van der Waals interactions in the MN_x structures. To visualize these results in real space, we have plotted three-dimensional representations with an isosurface value of s = 0.5 a.u. using the VMD software68 in Fig. S6 in the supplementary material. In this figure, green regions correspond to the van der Waals interactions between poly-N42− chains in MN_x. In addition, we can see some blue regions around metal atoms, which correspond to the attractive interaction between metal atoms and poly-N42− chains.

Nitrides containing poly-N42− chains have been synthesized for potential applications as HEDMs.34–37 This inspires us to study the detonation performance of the MN_x structures predicted in this work, which have more poly-N42− chains than the previously proposed structures. We assume that MN_x will decompose into MN and N₂ gas under ambient conditions: MNx(s)→MN(s)+12(x−1)N2(g). The estimated detonation properties of MN_x are listed in Table II. In particular, the energy densities of MN_x exceed that of TNT (4.30 kJ/g), except for P2₁-GaN₆ and P-1-YN₆. We calculate their detonation pressures P and detonation velocities V using the Kamlet–Jacobs empirical equations22,69P=15.58ρ2NM0.5Eg0.5 and V=1.01(NM0.5Eg0.5)0.5(1+1.30ρ), where N is the number of moles of N₂ gas generated during the decomposition of 1 g of MN_x and M is the molar mass of N₂ (28 g/mol). As can be seen from Table II, MN_x have excellent detonation pressures and detonation velocities that are higher than those of the explosives TNT and HMX.

By applying a structure search method accelerated by machine learning, we have successfully predicted four metal nitrides constructed from poly-N42− chains: P2₁-AlN₆, P2₁-GaN₆, P-1-YN₆, and P4/mnc-TiN₈. Based on first-principles calculations, we have found that the pressures required for synthesis of P2₁-AlN₆, P2₁-GaN₆, and P4/mnc-TiN₈ are around 40 GPa, and that required for P-1-YN₆ is near 20 GPa. Phonon calculations suggest that P4/mnc-TiN₈ is mechanically stable at ambient pressure, while the other three structures are mechanically stable at high pressures. Furthermore, the dynamical stabilities of MN_x at nonzero temperature have been shown by molecular dynamics simulations. Two types of poly-N42− chains are found: in type A, all the N atoms are sp²-hybridized and form delocalized π bonds distributed along the entire chain; in type B, half of the nitrogen atoms are sp³-hybridized, while the other half adopt sp² hybridization, and so π bonding can only occur for some of the N–N bonds in the poly-N42− chain. The results for DOS and pCOHP also confirm the MO analysis. The metallicity of MN_x stems from the delocalized π bonds. RDG calculations reveal that the effect of NCIs also helps to stabilize these MN_x structures. Most importantly, we estimate that the newly predicted MN_x have good detonation properties with potential applications as HEDMs. Furthermore, our predictions here can provide some guidance for future experiments.

Acknowledgment. J.S. gratefully acknowledges financial support from the National Natural Science Foundation of China (Grant Nos. 12125404, 11974162, and 11834006), and the Fundamental Research Funds for the Central Universities. K.X. acknowledges support from the National Natural Science Foundation of China under Grant No. 12004185, the Natural Science Foundation of the Jiangsu Higher Education Institutions of China under Grant No. 20KJB140016, and the Scientific Research Start-up Funds of Nanjing Forestry University (No. 163101110), and financial support from a Project funded by the China Postdoctoral Science Foundation (Grant No. 2019M651767). The calculations were carried out using supercomputers at the High Performance Computing Center of the Collaborative Innovation Center of Advanced Microstructures, the High-Performance Supercomputing Center of Nanjing University, and the High-Performance Computing Facility of Nanjing Forestry University.