Polymeric nitrogen containing single and double nitrogen bonds releases a large amount of energy during its decomposition and is thus a candidate green high-energy-density material (HEDM). However, the synthesis of stable HEDMs containing breakable single or double nitrogen bonds poses a real challenge. To date, there has been experimental confirmation of four polymeric nitrogen structures: cg-N (synthesized at 2000 K and 110 GPa),1 LP-N (3000 K and 150 GPa),2 HLP-N (3300 K and 244 GPa),3 and bp-N (2200 K and 146 GPa).4 However, these structures are synthesized under high temperatures and pressures, which limits their application as HEDMs. Therefore, it is an urgent, although formidable, challenge to explore novel polymeric nitrogen structures that can be synthesized under mild conditions.

Numerous theoretical and experimental results have indicated that the addition of metallic elements (M) to the nitrogen materials could provide new compound forms of polymeric nitrogen (MN_x) at lower pressures. Compared with pure polymeric nitrogen structures, these may have higher stability and milder synthesis conditions. Charge transfer from metal to polynitrogen may play important roles in reducing the pressure required for synthesis and improving structural stability. Theoretical calculations have predicted a variety of polynitride compounds MN_x (M = Li, Be, Na, Mg, Al, K, Ca, Cs, Rb, Ba, Se, He, or Xe, with x = 3, 5, 6, 10) featuring a variety of N₅ or N₆ polymeric nitrogen chains5–30 and some complicated nitrogen clusters [fused N₁₈ rings in KN₈ (Ref. 8) and N₁₀ rings in BeN₄ (Ref. 29), etc.]. The transition metals (Sc,31,32 Mo,33 Tc,34 Re,35 Fe,36,37 Ru,38 Os,39 Rh,40 Ir,39 Pd,41 Pt,41 Cu,42 Ag,43 Au,44 Zn,45 and Hg46) can also form nitrogen-containing structures. Some of these compounds are not potential HEDMs, but they do prove that introducing metallic elements to nitrogen seems to be a feasible way to enhance metastability. Experimentally, some metals or their compounds have been used as raw materials to synthesize polynitrogen compounds at relatively low pressures. CsN₅ has been successfully synthesized by compressing CsN₃ with N₂ near 60 GPa,47 and solid LiN₅ solid has been obtained by compression and laser heating of Li embedded in molecular N₂ at 45 GPa.48 Fe₃N₂, FeN₂, and FeN₄ have been synthesized from Fe and N₂ under different pressures,49 and elemental Cu reacts with N₂ to form CuN₂ at 50 GPa.50 Very recently, a beryllium tetranitride tr-BeN₄ in the triclinic phase has been synthesized at 85 GPa,51 and it has similar structural units as the previously predicted P-1 BeN₄ (infinite N chains coordinating Be atoms), proving that crystal structure prediction (CSP) searches can provide credible results.

The alkaline earth metals Be, Mg, Ca, Sr, and Ba are adjacent to the alkali metals in the Periodic Table, with similar properties but with a wider scope for chemical reactions and compound formation because of their additional valence electron. Among the alkaline earths, Be has the lowest atomic weight, and therefore its compounds with nitrogen have a high mass density of N and are good candidate HEDMs. A layered material, BeN₄, containing corrugated N₁₀ rings has recently been synthesized.51 Such a single-bonded polymeric nitrogen network has been predicted theoretically.29,52 Therefore, we have used an unbiased structural search with a particle-swarm optimization (PSO) algorithm53,54 to predict new beryllium nitrides with 1:4 compositions and have thereby found a new BeN₄ compound that may be stable at ambient pressure. In this compound, a planar N₄²⁻ ring is surrounded by four Be atoms, and the six π electrons of this ring conform to the Hückel [4n + 2] rule (with n = 1), indicating aromaticity of BeN_4. The large mass ratio of nitrogen in this compound (86.15%) leads to a high energy density (3.924 kJ/g).

We used the Crystal Structure Analysis by Particle Swarm Optimization (CALYPSO) code53,54 to search the structures of BeN₄ at different pressures, namely, 0, 10, 20, 50, and 100 GPa, at 0 K. Geometrical optimizations, total energy calculations, and electronic structure calculations were carried out in the framework of density functional theory (DFT) using the generalized gradient approximation with the Perdew–Burke–Ernzerhof functional (GGA-PBE),55 with a cutoff energy of 500 eV and a 0.03 Å⁻¹k-point grid56 in the Vienna Ab initio Simulation Package (VASP).57 The convergence threshold was used to optimize the atomic positions and lattice constants until the forces on each atom were less than 0.01 eV/Å and the energy change was less than 1.0 × 10⁻⁶ eV/atom. Projector augmented wave (PAW)58 pseudopotentials were adopted for Be and N, with 1s²2s² and 2s²2p³ valence states, respectively. Phonon dispersion relations were calculated using a supercell method with the finite displacement approach59 through the Phonopy code.60Ab initio molecular dynamics (AIMD) simulations were performed within the NVE ensemble61,62 to study thermal stability. The strain–stress methodology was used to calculate the elastic constants. The bulk and shear modulus, Young’s modulus E, and Poisson’s ratio v were obtained from the Voigt–Reuss–Hill approximations.63 The detonation velocity and pressure were calculated through the Kamlet–Jacobs semiempirical equations V_d = 1.01(NM^0.5E_d^0.5)^0.5(1 + 1.30ρ) and P_d = 15.58ρ²NM^0.5E_d,64,65 where N is the number of moles of gas per gram, M is the molar mass of the N₂ gas (28 g/mol), E_d is the energy density, and ρ is the mass density of BeN₄. The crystal orbital Hamilton population (COHP) was calculated to characterize the bonding properties using the LOBSTER program.66 Charge transfer was based on a Bader analysis,56,67 and the chemical bonding mode was calculated by the Solid State Adaptive Natural Density Partitioning (SSAdNDP) method.68

In a previous study, Zhang et al.29 used an unbiased structure search method to explore the structural evolutionary behaviors of beryllium polynitride compounds below 200 GPa. Interestingly, the only stable N-rich composition was found to be BeN₄, with the compounds P1-BeN₄ and P2₁/c-BeN₄ being theoretically predicted to exist at high pressure. In our work, also using an unbiased structure search method with CALYPSO, we have predicted a new BeN₄ in the tetragonal phase with the high group symmetry of P4/nmm. The calculated enthalpies of these three structures at pressures from 0 to 130 GPa at T = 0 K are shown in Fig. 1.

BeN₄ undergoes structural phase transitions first from P4/nmm to P1 at a pressure of 1.36 GPa and then to P2₁/c at 117.5 GPa. The phase transition pressure from P1 to P2₁/c phase that we calculated is similar to the reported value (115.1 GPa) in Fig. 1, which indicates that our method is appropriate. The predicted crystal structure of P4/nmm and its corresponding lattice parameters and atomic coordinates are shown in Fig. 2 and Table I, respectively. According to Fig. 2, there is a square planar N₄²⁻ ring in P4/nmm-BeN₄, where each inner N atom is shared by one Be atom and two adjacent N atoms. Within the cyclo-N₄, the length of the N–N bond is 1.363 Å and the angles of the N–N–N bonds are all around 90° at ambient pressure. Be atoms are situated at the center or vertex positions, with each Be atom binding to four adjacent N atoms with equal bond lengths of 1.713 Å. A variety of N₄ clusters have been predicted theoretically,24,30,69–71 including N₄⁴⁻, N₄²⁺, N₄²⁻, and N₄⁺. However, the novel square planar N₄²⁻ ring predicted in the present work has not previously been reported, although a square planar N₄²⁻ ring was found in bipyramidal Li₂N₄ by van Zandwijk et al.,72 and a series of alkali metal compounds M₂N₄ and alkaline earth compounds MN₄ with D_4h bipyramidal, D_2h planar, and C_2v planar structures were also found to contain such an N₄²⁻ ring.70,71,73 Here, we shall now focus on the square planar N₄²⁻ ring in P4/nmm-BeN₄.

Because there are no imaginary frequencies in the entire Brillouin zone (Fig. 3), the calculated phonon dispersions show that P4/nmm-BeN₄ is dynamically stable at ambient pressure. In addition, the phonon density of states shows that the low-frequency vibrational mode can be mostly attributed to the strong coupling vibration between Be and N atoms, while the N–N stretching mode is responsible for the high-frequency vibrational mode.

We further explore the thermal stability of BeN₄ at ambient pressure and high temperature by performing AIMD simulations on a 3 × 2 × 3 supercell containing 180 atoms. The images of the geometrical structure in Figs. 4(a) and 4(b) clearly reveal that the square planar N₄²⁻ ring maintains its structural integrity without any visible distortion after heating at 300 or 500 K for 10 ps with a time step of 1 fs. The radial distribution function (RDF) of N–N separations in P4/nmm-BeN₄ is shown in Fig. 4(c). The first sharp RDF peak is the average distance between the selected atom and its nearest surrounding atoms, here corresponding to the N–N bonds of BeN₄. This indicates that the N–N bond lengths remain at 1.369 Å as the temperature rises from zero to 500 K, which proves the stability of P4/nmm-BeN₄ at high temperatures. Thus, BeN₄ has a moderate cohesive energy and good thermal stability.

We also calculated the elastic constants of the P4/nmm-BeN₄ structure because these are helpful for determining the mechanical stability and hardness of a material. The elastic constants C_ij, bulk modulus B, shear modulus G, the Young’s modulus E and Poisson’s ratio v of P4/nmm-BeN₄ are presented in Table II. They are related by E = 9BG/(3B + G) and v = (3B − 2G)/[2(3B + G)].74 We calculated the bulk modulus of the β-Be₃N₂ hexagonal structure at ambient pressure to be B = 237.18 GPa, which is similar to a previously published value for this compound (236 GPa),75 which indicates that the calculation method we have adopted is appropriate for the Be–N system. The Pugh ratio k = B/G is used to estimate ductile or brittle characteristics, with k < 1.75 corresponding to brittle behavior and k ≥ 1.75 to ductile behavior. The value of k for the P4/nmm-BeN₄ structure is 1.02, indicating that this compound is brittle. For a tetragonal structure to be stable, the C_ij should satisfy the Born–Huang stability criteria C₁₁ > 0, C₃₃ > 0, C₄₄ > 0, C₆₆ > 0, C₁₁ − C₁₂ > 0, C₁₁ + C₃₃ − 2C₁₃ > 0, and 2(C₁₁ + C₁₂) + C₃₃ + 4C₁₃ > 0. The elastic constants of P4/nmm-BeN₄ satisfy these criteria, indicating that it should be mechanically stable under ambient conditions.

Energy density E_d is important for energy storage capability, while detonation velocity V_d and pressure P_d are known to be significant indices of detonation performance. The dissociation energy of P4/nmm-BeN₄ is calculated for the following decomposition path under ambient conditions: 3BeN₄(s) → Be₃N₂(s) + 5N₂(g). We use E_d = (96.4853 × BeN₄ molecule energy released)/(molar mass of BeN₄) to calculate the energy density of P4/nmm-BeN₄, the calculated E_d is about 3.92 kJ/g with reference to the lowest-energy phase at ambient pressure29,76 (the phase transition in Fig. S1 in the supplementary material). Table III compares the detonation properties of P4/nmm-BeN₄ with the experimental values for TNT and HMX. Although P4/nmm-BeN₄ shows superiority only with regard to its detonation velocity being higher than that of TNT, its detonation performance does indicate that it has potential application as an explosive.

To illustrate the electronic properties of BeN₄, its band structure, projected density of states (PDOS), and crystal orbital Hamilton populations (COHPs) under ambient conditions are presented in Fig. 5. There are no bands crossing the Fermi level in the band structure in Fig. 5(a), which indicates that the predicted BeN₄ crystal is a semiconductor with a bandgap of 3.67 eV. The obvious overlap between Be 2p and N 2p in the PDOS [Fig. 5(b)] reveals that the orbitals are strongly hybridized. Therefore, the bonding between Be and N is partially covalent. To characterize the orbital interaction between two N atoms, we performed a COHP analysis to measure the overlap strength between two N atoms. Negative values correspond to bonding states and positive values to antibonding states in Fig. 5(c), where it can be seen that the bonding states between N and N in BeN₄ are completely occupied and the antibonding states partially occupied, resulting in strong covalent bonding. The integrated COHPs (ICOHPs) calculated from the wave function’s energy-weighted population between two atomic orbitals represent the bonding strength. The ICOHP values of N–N and Be–N bonds are −12.21 and −1.63 eV, respectively. Therefore, the N–N interaction is substantially stronger than the Be–N interaction, which indicates that the covalent bonding between Be and N is weaker, which is consistent with the results of our previous hybridization investigation. Moreover, the Bader charges of Be and N atoms shown in Table S1 in the supplementary material reveal that the Be atoms combine with the N atoms primarily through electrostatic forces. Thus, the Be–N bonding has both covalent and ionic character.

We analyzed the chemical bonding patterns in the P4/nmm-BeN₄ crystal using SSAdNDP software. In a primitive cell, Be atoms connect four neighboring N₄ rings via four two-center, two-electron (2c-2e) Be–N σ bonds [Fig. 6(a)], and four N atoms form a ring through four 2c-2e N–N σ bonds [Fig. 6(b)]. There are two types of π orbital [Figs. 6(c) and 6(d)]: one links two adjacent N atoms through two 2c-2e N–N π bonds, while the other links four N atoms via one four-center, two-electron (4c-2e) N–N large π bond in a plane. Together, the σ and π bonds require 2 × 11 = 22 electrons, which is equal to the number of outer valence electrons of BeN₄. It can be seen that the localized 2c-2e σ bonds, 2c-2e π bonds, and delocalized 4c-2e π bonds have an influence on the high structural stability of BeN₄. The total of six π electrons is the most striking feature, since it satisfies the Hückel [4n + 2] rule (with n = 1), indicating that BeN₄ is an aromatic structure.

The band structure, PDOS, COHP calculations, and chemical bonding pattern analysis illustrate the covalent bonding between the N atoms. The Be atoms contribute their outer-shell electrons to the N₄ cluster and use the outer-shell 2p orbital to accommodate the lone-pair electrons of N₄²⁻, with the formation of covalent and ionic bonds to stabilize the crystal. This explains why Be can stabilize the square N₄ ring.

Be₃N₂ with two known phases (α80,81 and β82) has been experimentally synthesized. We have calculated the enthalpy of formation of the predicted P4/nmm-BeN₄ relative to Be₃N₂ and N₂ (the phase transition in Fig. S1 in the supplementary material) to find a possible high-pressure synthetic path. Based on the results in Fig. 7, the predicted BeN₄ can be synthesized through mixing Be₃N₂ and N₂ together at pressures above 31.6 GPa, which indicates that square planar N₄²⁻ can be synthesized beyond 31.6 GPa. It should be noted that both the P2₁/c and P-1 phases have lower energies than the P4/nmm phase above 31.6 GPa (Fig. 1), which makes the synthesis of P4/nmm-BeN₄ a difficult task. Fortunately, a number of experimental results achieved with high-pressure techniques can provide inspiration for the synthesis of P4/nmm-BeN₄. The recently synthesized bp-N is not the phase with the lowest energy at the corresponding pressure.4,83 Similarly, silicon and ice can form different phases along different pressure–temperature paths.84,85 We anticipate that experimental investigations will also find the precise pressure and temperature conditions for the synthesis of P4/nmm-BeN₄.

Through structural searches and density functional calculations under ambient conditions, we have found a P4/nmm-BeN₄ structure with high kinetic and thermodynamic stability. The phonon dispersion and electronic properties of this structure indicate the presence of covalent bonding between two N atoms, with covalent and ionic bonding between Be²⁺ and square planar N₄²⁻ enhancing the stability of the semiconductor BeN₄ structure. The total of six π electrons satisfies the Hückel [4n + 2] rule (with n = 1), indicating that BeN₄ has an aromatic structure. Furthermore, the energy density of P4/nmm-BeN₄ is 3.92 kJ/g, making it a candidate high-energy-density material. The calculated formation of enthalpy of BeN₄ indicates that the square planar N₄²⁻ ring might be synthesized through mixing Be₃N₂ and N₂ at 31.6 GPa. The results of this investigation of the role of novel chemical bonds in nitrogen chemistry should encourage experimental efforts to synthesize these promising high-energy materials.

Acknowledgment. This work was supported by the National Natural Science Foundation of China under Grant Nos. 11974154, 11674144, and 11604133, and the Natural Science Foundation of Shandong Province under Grant Nos. ZR2018MA038, 2019GGX103023, Z2018S008, and 2019KJJ019.