Crystalline Structures and Energetic Properties of Lithium Pentazolate under Ambient Conditions

Recently, it has been reported that high-pressure synthesized lithium pentazolates could be quenched down to ambient conditions. However, the crystalline structures of LiN5 under ambient conditions are still ambiguous. In this work, the structures of LiN5 compound were directly explored at atmospheric pressure by using a new constrain structure search method. By using this method, three new allotropes were confirmed, and they show lower energy than the previous reported LiN5 phases. Both their thermodynamic and dynamic stability were confirmed through formation enthalpies, phonon spectrum, and ab initio molecular dynamics simulations under ambient conditions. Moreover, these three allotropes show similar formation enthalpies and properties, which suggests that it is hard to obtain a single LiN5 phase, which is well consistent with the experimental phenomenon. Furthermore, because of their low formation energy, all of them possess low energy density when they directly decompose to Li3N and nitrogen (0.52 kJ/g). Instead, the decomposed energy could be further improved to 3.78 kJ/g when they decompose under an oxygen-rich environment.


INTRODUCTION
Pentazolate anion (cyclo- N 5 − ) is receiving ever-increasing significant interest because of its potential application in highenergy density materials. 1−9 Especially, since 2016, several breakthroughs have been witnessed in this field, including some metal-free pentazolate hydrates [such as (N 5 ) 6 6 (N 5 ) 2 ]· 4H 2 O} synthesized in experimentation under ambient conditions. 10−21 However, the nonenergetic group (such as water and halogen) lowers the nitrogen concentration as well as energy density of these compounds. 18,22 One important way to obtain nonenergetic group metal pentazolates (such as CsN 5 , NaN 5 ) is by using high pressure and laser-heating, but few of them could be quenched down to ambient. In 2015, Peng et al. and Shen et al. independently predicted that LiN 5 could be synthesized under high pressure (15−100 GPa) through theoretical prediction, 23,24 and Laniel et al. successfully achieved this goal by the experiment method at 2018. 25 Moreover, the synthesized compounds could remain stable when releasing pressure to ambient conditions. These reports are important for synthesizing and applying polynitrogen materials under atmospheric pressure. However, Laniel et al. could not obtain a single clean diffraction pattern containing solely the diffraction lines despite many attempts, and the crystalline structure of the LiN 5 compound is still unclear. 25 There are series of excellent reports that worked on predicting metal pentazolates, which usually can become thermodynami-cally stable under high pressure, such as CuN 5 , LiN 5 , NaN 5 , ZnN 10 , and AlN 15 . 23,24,26−31 However, these predictions are hard to be directly applied to the ambient conditions, which require to directly search metastable structures. Notably, it is still a big challenge to directly find the most stable microstructures of metal pentazolate compounds by the pristine structure search algorithm because N 2 molecules have much more selective advantages than the N 5 ring in energy. Also, the structures were obtained from high-pressure search, and there is no guarantee that they have the lowest energy at atmospheric pressure and very often they do not. 24,32 Recently, we have overcome this difficulty by proposing a constrained crystal structure search method to directly predict the most stable MgN 10 , AlN 15 , and CuN 5 crystalline structure and property under ambient conditions. 33 In this work, to answer the above questions, the constrained crystal structure search method was used to explore the crystalline structures between 0 and 15 GPa. Here, three possible new allotropes of LiN 5 under the ambient condition were reported, and the structure phases between 0 and 15 GPa were confirmed in this work. Moreover, the stabilities of these structure phases were evaluated by formation enthalpies, phonon spectrum, and ab initio molecular dynamics (AIMD) simulations. These structures have a great possibility to coexist at a certain low pressure because they have very close formation enthalpy and similar properties, which could be used to explain that why Laniel et al. could not obtain solely the diffraction lines of the LiN 5 compounds. Furthermore, these compounds show low energy density (about 0.53 kJ/g) when directly decompose to Li 3 N and N 2 . An alternative approach is to make LiN 5 decomposing under an oxygen-rich environment, and the energy density would improve to 3.78 kJ/g.

CALCULATIONAL METHODS
The structure predictions were through a global minimization of potential energy surfaces based on the particle swarm optimization methodology as implemented in the CALYPSO code. 36−38 The stoichiometric (LiN 5 ) x (x = 2−6) up to 3000 structures was directly searched in our simulations under 1 atm. The cyclo- N 5 − was constrained as a group unit to produce initial structures, and all the structures were optimized by the local density functional theory (DFT) soft package until converged. Also, then, structures were ranked according to their enthalpy, and the top 50 structures were further optimized with finer convergence criteria and used for further stability analysis.
The local structural relaxations and electronic properties calculations were performed in the framework of DFT as implemented in the Vienna Ab initio Simulation Package (VASP). 39,40 The one-electron wave function was expanded using a plane-wave basis set with an energy cutoff of 750 eV. The projector augmented wave method 41 was used to describe the interactions between the ion and electron. Also, electron correlation was treated with the generalized gradient approximation with the Perdew−Burke−Ernzerhof (PBE) functional, 42 in which 1s 2 2s 1 and 2s 2 2p 3 are treated as valence electrons for Li and N atoms, respectively. The atomic positions and lattice constants were optimized using the conjugate gradient scheme until the residual Hellmann−Feynman forces on each atom were less than 0.01 eV/Å, and the self-consistent field calculations were stopped until energy change was smaller than 1 × 10 −5 eV/atom. Brillouin zone (BZ) integrations are performed using a Monkhorst−Pack 43 k-point mesh with a resolution of 2π × 0.03 Å −1 . The phonon spectrum was calculated by the density functional perturbation theory method with a supercell containing about 100 atoms for different phases, as implemented in Phonopy. 44 The AIMD simulations were implemented using an NVT ensemble with a large supercell containing about 100 atoms for different phases to evaluate their thermal stability. The Raman spectrum was calculated by CASTEP code, and the Raman peak (1331 cm −1 ) of diamond was used for peak position corrections. 45 The convergence test results are shown in Figure S1, which prove that our computational parameters are reliable. The formation enthalpy (H f ) was calculated by using the following formula (refers to the reaction of 2·Li + 5·N 2 → 2LiN 5 ) in which enthalpy H is obtained for the thermodynamically stable structures of certain compositions at each given pressure. For pure N 2 , Pa3-phase, P42/mnm-phase, and P2 1 /c-phase nitrogen crystal were considered. 46 For pure Li, the bcc-phase, fcc-phase, cI16-phase, and Aba2-40-phase were considered. 47 The decomposed energy (E d ) of LiN 5 compounds was evaluated by directly decomposing with the reaction of 3·LiN 5 → Li 3 N + 7· N 2 by using the following formula (2) in which energy H is obtained for the thermodynamically stable structures of certain compositions at 1 atm, and M LiN 5 is the relative atomic mass of LiN 5 . Also, under an oxygen-rich environment, with the reaction of 4·LiN 5 + O 2 → 2·Li 2 O + 10· N 2 was used to calculate the E d by using the following formula (3) All the calculations were assisted by qvasp. 48 3. RESULTS AND DISCUSSION 3.1. Formation Enthalpy and Energy Density. After thousands of times calculations and data screening, we get three new LiN 5 allotropes, whose formation enthalpies are both negative at 1 atm. Their space groups are P6 1 22, P2 1 2 1 2 1 , and Pnan, so they were named P6 1 22-LiN 5 , P2 1 2 1 2 1 -LiN 5 , and Pnan-LiN 5 , respectively. As shown in Figure 1a, we listed the formation enthalpies of these seven LiN 5 allotropes, including three new phases and four have already reported phases (obtained from high pressure) (named P2 1 -LiN 5 , P2 1 /m-LiN 5 P2 1 /c-LiN 5 , and C2c-LiN 5 ) 23,24 in the pressure range of 0−15 GPa. Notably, the four allotropes, including P6 1 22-LiN 5 , P2 1 2 1 2 1 -LiN 5 , Pnan-LiN 5 , and P2 1 /c-LiN 5 , are more stable than the others because they have negative formation enthalpy at 1 atm (shown in Table 1). If ignoring the kinetic barrier, they show great possibility to be formed when their high-pressuresynthesized compounds were released pressure to ambient conditions. In particular, the formation enthalpy curve of P6 1 22-LiN 5 decreases very slowly when compared to other allotropes, which indicate that P6 1 22-LiN 5 is harder to directly form under high pressure. Furthermore, the formation enthalpy curve of P2 1 2 1 2 1 -LiN 5 and P2 1 /c-LiN 5 is very close at 0−15 GPa (small than 0.07 eV/atom), and these two phases are with great possibility to coexist at low pressure that might be used to explain why Laniel et al. could not obtain solely the diffraction lines of the LiN 5 compounds. 25 Also, the formation enthalpies at higher pressure of up to 100 GPa of these seven allotropes are listed in Figure S2, and as shown in Figure 1a, we add structural phase information at an unknown range of 0−15 GPa. Furthermore, as shown in Figure  1b, we confirmed that the P6 1 22-LiN 5 has lower formation enthalpy at 0−1 GPa, P2 1 2 1 2 1 -LiN 5 occupied at 2−14 GPa, P2 1 / c-LiN 5 occupied in 14−15 GPa, P2 1 /m-LiN 5 occupied in 15−62 GPa, and C2/c-LiN 5 at 62−100 GPa. Then, we list the decomposed energy of LiN 5 with and without oxidant, as shown in Table 1. We can directly see the remarkable contradiction between high energy density and high stability, and the P6 1 22-LiN 5 shows the lowest energy density because of low formation enthalpy. All the LiN 5 allotropes show low energy density (0.52−1.38 kJ/g) when directly decomposed under an oxygen-free environment. Also, the energy density improved to 3.68−4.55 kJ/g when directly decomposed under an oxygenrich environment, which is higher than TNT. Considering the explosion pressure obtained from large amounts of gas N 2 molecules released from an oxidizing reaction, the compounds of LiN 5 are still great potential application in propellant and explosives.
3.2. Structure Features. As discussed above, here, we focus on the structures of LiN 5 compounds with negative formation enthalpy, which include P6 1 22-LiN 5 , P2 1 2 1 2 1 -LiN 5 , Pnan-LiN 5 , and P2 1 /c-LiN 5 , because they have great possibility to maintain stability at ambient conditions. The lattice constant and Wyckoff positions of them are shown in Table S1. As shown in Figure 2a,  5 . Also, there are three type Li−N bonds, including the length of 2.31, 2.10, and 2.14 Å (Figure 2a), which is a bit larger than the average bond length (2.11 Å) of Li−N in P6 3 /mmc-Li 3 N.
For P2 1 2 1 2 1 -LiN 5 , P2 1 /c-LiN 5 , and Pnan-LiN 5 , they are all composed of the tetrahedral Li configuration unit, and each Li forms four bonds with the η 4 -N 5 ring. Also, they just show different space groups. As shown in Figure 2b

ACS Omega
http://pubs.acs.org/journal/acsodf Article result in very close formation enthalpy between these three phases at ambient conditions. Suffering from less coordination number, these three phases have a bit higher formation enthalpy (0.072 eV/atom) than P6 1 22-LiN 5 , and these three phases show more homogeneous bonding features, which increases their possibility of coexisting at ambient conditions. 3.3. Stability Evaluations. Furthermore, the stabilities of these allotropes were evaluated by dynamic and thermodynamic  Furthermore, we used AIMD simulations to evaluate the thermodynamic stability of these allotropes. As shown in Figure  4a−d, the total energies of these allotropes both fluctuate around a certain energy level, the fluctuations are quite small, and the images of the geometry structure at the end of simulations clearly reveal that the structural skeleton did not suffer a large shape change up to 600 K. When temperature increases to 800 K, as shown in Figure S3a−d, the total energies experience violent vibration around the reference energy line, and the structures at the end of simulations clearly revealed that the cyclo-N 5 − maintains its structural skeleton, except for clear change for Li tetrahedron. When temperature continues to increase to 1000 K, cyclo-N 5 − decomposed to N 2 molecules. The simulation results were consistent with our previous reports 33 that the N−N bonds are harder to be decomposed than the metal−N bonds, suggest that it is much easier to isolate cyclo- The charge information on each atom was further calculated by Bader's quantum theory of atoms in molecules 50,51 and is shown in Figure 5e−h, and we can see that the Li atoms play a role of the electron-donating unit, and each Li atom loses about 0.86 electrons to cyclo-N 5 − and shows +1 valence state. The charge distribution on cyclo-N 5 − is more homogeneous (0.10− 0.20) than many known pentazolates ( Figure S4), which reveals good aromaticity in these compounds. Especially in P6 1 22-LiN 5 , each Li atom has five coordination number, every N atom in cyclo-N 5 − is coordinated with the Li atom, and the electrons on the N atom in cyclo- N 5 − are roughly equal, which caused better aromatic property and lower formation enthalpy than other allotropes.
We employ crystal orbital Hamilton population (COHP) to analyze the bond strength quantitively. COHP is an energyweighted DOS, which describes bonding and antibonding state arrangements in the band energy scale. Moreover, −ICOHP was used to estimate the overlap strength of N−N bonds and Li−N bonds, which was calculated by using the LOBSTER program. 52 The larger positive value means stronger bonding between two atoms. As shown in Figure 5i  The calculated electronic band structure of LiN 5 allotropes is shown in Figure 6a−d. The results reveal that these two LiN 5 compounds are both semiconductors with a large band gap (both than 5 eV) at the PBE level at 1 atm. Considering density functional calculations usually lead to a considerable underestimation of the energy gap, and the actual band gaps are expected to be larger than the calculated results, thus the LiN 5 compounds would present the transparent crystal at ambient conditions because large optical band gap will induce weak visible light absorption. 34 The band gap as a function of pressure is shown in Figure S5 from the range of 0−15 and 0−100 GPa. The band gap of all LiN 5 allotropes decreases with the increase of pressure, and they still show a big band gap at high pressure except P6 1 22-LiN 5 , which is consistent with the experimental phenomenon that microphotographs appear translucent at 52.0 GPa. 25 The Raman spectrum of experimental data and all the LiN 5 allotropes under ambient conditions is shown in Figure 7, and more detailed low-intensity Raman peaks are shown in Figure  S6. The calculated results reveal that the breathing and stretching cycle-N 5 − Raman modes could still be measured at peaks around 1200, 1130, and 1035 cm −1 in all LiN 5 allotropes. The faint low-frequency peak (364 cm −1 ) is observed in P2 1 2 1 2 1 -LiN 5 , P2 1 /c-LiN 5 , and Pnan-LiN 5 (both 374 cm −1 ), and these results are well consistent with the experimental. 25 Moreover, many newly observed low-intensity peaks are marked in Figure 7, and except P6 1 22-LiN 5 , all other allotropes are possible to exist or coexist in experimental synthetic samples, which could be used to make clear that why Laniel et al. cannot get solely the diffraction lines of the LiN 5 compounds despite many attempts.

CONCLUSIONS
In summary, here, we reported three new LiN 5 compounds with negative formation enthalpy and supplemented the structure phases between 0 and 15 GPa. The stability of them is both confirmed by the phonon simulations, AIMD simulations, and bonding analysis. The calculated results reveal that the N−N bonds are harder to decompose than Li−N bonds, and the cyclo-N 5 − can be maintained at least to 600 K, which is enough to be as

■ ACKNOWLEDGMENTS
We sincerely thank Prof Maosheng Miao for his helps and discussions. This work was supported by the National Natural Science Foundation of China (grant nos. 21905159, 11974207 and 11974208), Natural Science Foundation of Shandong Province (grant nos. ZR2019BA010, ZR2019MA054, and 2019KJJ020), and the Project of Introduction and Cultivation for Young Innovative Talents in Colleges and Universities of Shandong Province. Calculations were performed at the High-Performance Computing Center (HPCC) of Qufu Normal University.