Predicting high energy density materials using ab initio evolutionary algorithms: The results for N 5 AsF 6

N 5 AsF 6 is the ﬁrst successfully synthesized salt that has a polymeric nitrogen moeity (N + 5 ). Although 12 other N + 5 salts followed, with N 5 SbF 6 and N 5 Sb 2 F 11 being the most stable, the crystal structure of N 5 AsF 6 remains unknown. Currently, it is impossible to experimentally determine the structures of N 5 AsF 6 due to its marginal stability and explosive nature. Here, following an ab initio evolutionary prediction and using only the stoichiometry of N 5 AsF 6 as a starting point, we were able to reveal the crystal structure of this high energy density material (HEDM). The C 2 V symmetry of the N + 5 cation, as suggested from earlier investigations, is conﬁrmed to be the symmetry adopted by this polymeric nitrogen within the crystal. This result gave full conﬁdence in the validity of this crystal prediction approach. While stability of the N + 5 within the crystal is found to be driven by electronic considerations, the marginal stability of this HEDM is found to be related to a partial softening of its phonon modes. conducted using density functional perturbation theory 37 with PBEsol. The electronic Brillouin zone integration as well as phonon calculations were performed on a Monkhorst-Pack 12 × 12 × 12 k-point meshes and a 2 × 2 × 2 mesh of phonon vectors for the dynamical matrix computation. The coordinates of the two molecular units making the crystal were fed to GamessUS 38 linked with nbo.6 39 to carry out the natural bonding orbital (NBO) theory analysis without boundary conditions and without any further geometry optimization. Partial charges were then extracted from this analysis. The true bonding within N + 5 was also extracted from the bond order obtained with this NBO analysis. The NBO investigation was performed at the PBE level of theory with PBEsol and B3LYP/6-31++G. Xcrysden 40 , VESTA 41 and Jmol 42 were used for the visualization of the results and the plots of this manuscript.


Introduction
Polymeric nitrogen (PN), a form of singly or doubly bonded nitrogen atoms as clusters or extended solid networks, has been the subject of plethora of theoretical and experimental investigations. The interest in these materials originate from their energy storage capabilities. With bonding energies at: 160 kJ.mole −1 for nitrogen single bond N-N, 418 kJ.mole −1 for nitrogen double bond N=N and 954 kJ.mole −1 for nitrogen triple bond N≡N 1 , an energy release of 794 kJ.mole −1 could be harnessed from the single to the molecular triple bond transformation. This is equivalent to tenfold higher detonation pressures compared to the high energy explosive HMX 2, 3 . The pioneering work of McMahan and Lesar 4 in 1985 stands among the first attempts aiming at investigating the existence of polymeric nitrogen structures. They predicted that the transition pressures from molecular to monoatomic nitrogen was less than 1 Mbar. This was well within reach of most high pressure anvil cells in use and the simple cubic structure formed showed the common distortions found in group V elements. In 1992, McMahan et al. 5 theoretically identified more extended nitrogen solids at relatively low pressures that are within reach of typical diamond anvil cells. Notable among these structures is the cubic gauche phase with three-fold coordinated nitrogen atoms adopting all-gauche dihedral angles. A 2004 experimental breakthrough by Eremets et al. 6 succeeded in synthesizing this structure by pressurizing molecular nitrogen above 110 GPa at temperatures above 2000 K by laser heating in a diamond anvil cell. Another polymeric nitrogen structure in the form of two N 8 isomers per unit cell held together by van der Walls forces was predicted by Hirshberg et al. 7 but has yet to be synthesized. The calculations showed that this polynitrogen form is stable even at ambient pressure and is more stable than the cubic gauche phase below 20 GPa. There is a strong debate whether the synthesis of polymeric nitrogen could be achieved beyond out-of-equilibrium processes, such as high pressure combined with high temperature and plasma methods 6,8 . To this date, extended solid forms of polymeric nitrogen were only obtained under extreme conditions and in most cases the structures were lost upon releasing pressure. To this date, extended solid forms of polymeric nitrogen were only obtained under extreme conditions and in most cases the structures were lost upon releasing pressure. Cluster polynitrogens on the other hand had substantial success with the pentazolium N + 5 cation stabilized as the N 5 AsF 6 salt in the work of Christe et al. 9 . This received a lot of attention in the scientific community as well as in many media outlets, such as the New York Time (see Supplementary Note-1) and also led to substantial sponsorship of polynitrogen chemistry 10 by the Defense Advanced Research Projects Agency (DARPA).
Marginal stability of the N 5 AsF 6 solid warrants explosive hazard possibilities that prohibits the use of X-ray diffraction to determine its crystal structure. Thus, no known crystal structure of this compound exists. Attempts to clarify the thermal stability of N + 5 salts can be found in the work of Yu 11 , where the relative stability of the N + 5 salts is attributed to the role played by the central atom of the counter anion (e.g. As, Sb. Al, B . . . etc) and their ligands (eg F, (OH) 4 F 2 , (CF 3 ) 4 . . . etc). To the best of our knowledge, the only known crystal structure of a N + 5 salt is that of N 5 Sb 2 F 11 , which was synthesized using the N 5 SbF 6 precursor 12 (The crystal arrangement of N 5 Sb 2 F 11 can be found in Supplementary Fig. 1). Fully understanding polymeric nitrogen salts such as N 5 AsF 6 will shed more light and guide the experimental search for the long-sought-after polymerization of nitrogen. For example, in light of the most recent synthesis of cyclic N − 5 in (N 5 ) 6 (H 3 O) 3 (NH 4 ) 4 Cl 13 , a N 10 synthesis from N + 5 and N − 5 precursors is worth investigating. In this paper, using an evolutionary ab initio search procedure, we were able to obtain the structure of N 5 AsF 6 . We therefore performed a thorough investigation into the overall stability of the crystal from a phonon perspective. We also assessed N + 5 stability within the crystal from an electronic perspective. Despite the tremendous success of evolutionary search techniques in the last 10 years in predicting superconductors [14][15][16][17] , magnetic materials 18,19 and novel compounds that defy conventional chemistry 20,21 , the work of Pakhnova et al. 22 stands out as the only predictive search fully dedicated to the study of energetic materials including most known explosives such as TNT, PETN, β -HMX, CL-20, TATB and their mixtures. Taking into account that clusters should retain chemical stability in the crystal as a requirement for the overall stability of the compound 7 , the crystal prediction methodologies will identify polynitrogens not only in the gas phase, but it will also help in identifying the best clusters and their embodiment with other elements across the periodic table. An evolutionary search, guided by quantum mechanics ab initio calculations, spans from hundreds to thousands of possible structures and survival of the fittest guides the search towards the true structure.

Structure search
The search for the crystal structure of N 5 AsF 6 was carried out using the evolutionary algorithm for structure predictions as developed in the code USPEX. A single unit with five nitrogens, one arsenic and six fluorine atoms and with no other chemical or structural knowledge set the basis for the calculations. In the course of a typical search, generations of structures were produced by the USPEX code and geometry optimization was conducted on each structure using density functional theory (DFT). As the first produced generations are more likely to be unphysical, geometry optimization often fails to converge. Additionally, nitrogen and fluorine in N 5 AsF 6 structures represent yet another challenge due to their comparable ionic size and electronegativity. Hence, computational compromises are needed to avoid too many failing structures leading to the halt of the evolutionary search. Another hurdle is the inherent marginal stability of the compound itself as reported in its synthesis and characterization procedures. Low temperature Raman and IR vibrational frequencies as well as DFT calculations hypothesized the C 2V symmetry of the N + 5 cation 9, 12 . In this work, the zero temperature, zero pressure evolutionary calculations succeeded in producing this C 2V polymerization of nitrogen in the best three structures as shown in Figs. 1-a, 1-b and 1-c; thus, representing a validation of the methodology used. In addition, the AsF6 − anion octahedral O h symmetry was also obtained in the best three predicted structures.
The nitrogen C 2V as well as the AsF − 6 configurations are lost in the fourth best structure that is shown in Fig. 1-d. In this structure, the system exhibits a decomposition into molecular N 2 , nitrogen trifluoride NF 3 and arsenic trifluoride AsF 3 . Moreover, the aforementioned fourth structure, where the N + 5 polymeric nitrogen is lost, is also far in its enthalpy from the first three structures. The energy shift of the fourth best structure is as high as 26.3 meV/atom from the best structure. The shift, however, is only 3.1 meV/atom and 5.1 meV/atom from the second and third best structures, respectively. The three best structures also set themselves as a separate subset of structures in the volume landscape with a volume of ∼16 3 /atom. This is among the lowest volumes found in comparison to all the structures obtained along the whole evolutionary search procedure with the USPEX algorithm. The rest of the unfavorable structures possess higher volumes than this special set. This suggests the stability of the N + 5 in this solid at relatively smaller volumes and confirms the well-known polymorphism of polymeric nitrogen structures like the singly bonded cubic gauche to form at higher pressures.
The best energetically favorable structure shown in Fig. 1-a (different orientations of the structure are shown in Supplementary Fig. 2 ), is found to be triclinic of space group H1 P1 (Crystallographic Information File (CIF) provided in Supplementary note-2) based on DFT calculations with PBEsol general gradient approximation (GGA). In fact, neither PBE GGA nor LDA approximations were capable of handling the evolutionary search. This is due to the computational difficulty involved in configuring an odd number of nitrogen atoms in the unit cell and nitrogen affinity to form molecular N 2 ; the dangling nitrogen left in the presence of the AsF − 6 anion thus becomes problematic in terms of stability and computational convergence. Using PBE GGA and LDA functionals, the evolutionary search often comes to end without achieving the task as the majority of the structures cannot satisfy the system constraints. Only in one instance was LDA able to finish the search but without achieving polymerization of nitrogen into N + 5 fragments (See Supplementary Fig. 3 for the best three structures obtained using the local approximation). While LDA lacks the delocalization required, PBE GGA overestimates the lattice parameters. PBEsol on the other hand is notorious for systematically lowering the lattice constants compared to other GGA functionals which improves equilibrium properties of densely packed solids and their surfaces 23 . With the volume shrinkage, polymerization of nitrogen becomes easier which explains PBEsol capability to obtain the desired N + 5 polymorph. The obtained structure was further relaxed in a variable cell geometry optimization with PBEsol, PBE and LDA. The results are summarized in Table 1 Similar search strategy was used for the N 5 SbF 6 and N 5 Sb 2 F 11 . Unfortunately the evolutionary search for N 5 Sb 2 F 11 was computationally too expensive and only a few structures satisfied the system constraints to progress towards better structures thus the halt of the algorithm. However, USPEX algorithm produced a similar triclinic structure for N 5 SbF 6 with N + 5 C 2V and AsF6 − O h symmetries preserved and with the following lattice constants: a = 18.894 , b = 17.262 , c = 16.917 , α = 85.73 • , β = 98.5760 • and γ= 89.233 • . The crystal arrangement of N 5 SbF 6 can be found in Supplementary Fig. 4.

Stability, vibrational spectroscopy and phonon calculations
There are 12 atoms in the unit cell for N 5 AsF 6 , the phonon spectrum thus contains 36 branches (3 acoustic and 33 optical). These 36 branches correspond to 36 normal vibrational modes at the center of the first Brillouin zone Γ. Computational spectroscopy to get vibrational modes using the linear response methodology is only performed at the Γ point. A tentative comparison with the observed IR signals regarding the N + 5 of reference 9 are summarized in Table. 2 and the total vibrational data is provided in Supplementary note-3. No negative frequencies are found upon imposing the acoustic sum rule. Also, excellent agreement can be found with the experimental Raman and infrared frequencies reported in the work of Christe 9 Table 2. Computational IR activity using the linear response approach at the PBEsol level of theory performed at Γ. Four N + 5 -related normal modes observed in the work of Christe et. al 9 are presented here for comparison. Frequency wavenumbers are in cm −1 and IR intensities in (D/A) 2 /amu capability in producing exact bond lengths for the N + 5 terminal. Mode 34 on the other hand is largely overestimated. This is mostly due to the strong anharmonicity within DFT especially for the ν 8 asymetric central stretch.   2 shows the phonon dispersion of N5AsF 6 based on the GGA approximation in the PBEsol formalism. We adopted the notations and crystallographic directions suggested in 24 for such a triclinic structure of type TRI1b with the k-path: X−Γ−Y|L-Γ−Z|N-Γ−M|R-Γ. Two partially unstable phonon modes as indicative from the negative frequencies developed can be observed along the high symmetry points of first the Brillouin zone. This softening is more severe towards Γ point with the largest modulus in the Z-Γ and N-Γ directions. It is worth noting that this phonon branch corresponds to a collective motion of the crystal units in the same direction. In light of this finding, a trivial solution would be a structural phase transition to a lower symmetry structure to obtain more stable phonons. However, the determined H1 P1 space group for the N 5 AsF 6 in this predictive search is already a low symmetry triclinic structure. Consequently, the marginal stability exhibited by N 5 AsF 6 can be explained by the existence of these somewhat soft phonon branches. The computed phonon of this branch that causes instability of the structure at 0 K corresponds to the librational motion of the overall crystal in the same direction that is about 45 degrees from the [010] direction as shown in Fig. 3. This displacement is responsible for the perturbed stability of this compound up to an unknown temperature that should be within the experimental conditions of the synthesis and characterization of N 5 AsF 6 . It is important to note that ab initio calculations are conducted at zero temperature and zero pressure conditions while the synthesis and the spectroscopic data of reference 9 were at ambient pressure and at temperatures of -130 • C for Raman and -196.8 • C for IR. Moreover, high energy density materials such as N 5 AsF 6 are highly explosive. Impact sensitivity and performance of such compounds are considered the key requirements towards their development. Upon impact, the first excitation takes place on the acoustic phonon at low wavenumbers, this causes energy to be transferred to higher energy phonon modes. For instance, bond breaking follows a stretching mode (about 1000 cm −1 and higher) of the compound units before detonation can occur 25,26 . Softening of low frequency phonons has been observed in other high energy density materials, such as thalium azide TlN 3 below 240 K 27 . Therefore, N 5 AsF 6 can be seen as a high explosive with higher sensitivity upon impact due to its thermally sensitive, soft acoustic phonons.

Electronic structure and NBO investigation
In Fig. 4, the electronic structure of N 5 AsF 6 is plotted along the high symmetry points of the first Brillouin zone. Similar to most energetic materials, the N 5 ASF 6 salt is an insulator. It has a large direct bandgap of 3.82 ev at the PBEsol level of theory. More importantly, Natural bonding orbital (NBO) analysis is more efficient in studying the intramolecular and intermolecular bonding that takes place at the molecular level. Hence its implementation in this investigation. N + 5 stability was long believed to be attributed to its structural and electronic resonance 9,12 . N + 5 and AsF − 6 as optimized in the crystal were fed to the GamessUS density functional theory package (2019 R1 version) in combination with the natural bonding orbital method as implemented in nbo.6 at the PBEsol level of theory. A summary of natural population analysis is shown in both Table 3 and Fig.5. This analysis provides the correct resonant N + 5 configuration (A total of six resonant structures are debated to explain N + 5 existence 9, 12 ). Partial natural charges on N + 5 as shown in Fig. 5 and Table 3 in our work are within 7-11% difference in comparison to the work of Fau and Bartlet 28 . This discrepancy is mainly attributed to the electronic effects imposed by the neighbouring AsF6 − as delivered from the solid in our investigation. NBO analysis also shows that N + 5 polymorphism is largely driven by electronic considerations. From bond order analysis, the cation is encapsulated into stability through the strong triply bonded nitrogens from both ends (N8-N9 and N11-N12 with about 2.3 bond order in Fig. 5) despite the fragile central single bonds (N9-N10 and N10-N11 with about 1.1 bond order in Fig. 5). For comparison, the well-known azide salts such as NaN 3 which represent a class of high energy density materials have the azide anion N − 3 stabilized through a similar encapsulation of the central atom electronically by resonant double and triple bonds within its linear chain.
Compared to the only known crystal structure among all N + 5 salts which is N 5 Sb 2 F 11 from reference 12 , excellent agreement was found for the intramolecular bond lengths and angles. A summary of this comprehensive comparison can be found in Table 3 with a focus on the N + 5 fragment only since the counter anion in N 5 Sb 2 F 11 is not the same and involves two SbF − 6 units sharing a bridging fluorine. Similar to end-on addition of N + 5 and N − 5 salts to obtain compounds with higher nitrogen content, earlier attempts to investigate addition of N + 5 and N − 3 salts was conducted by Fau and Bartlett 28 in the gas phase. Isolating a covalently bonded N 8 is found to be difficult because decomposition is energetically more favourable. However, we believe that taking the whole crystal of N + 5 into consideration is necessary since a large lattice energy within a solid-state structure was long believed to allow such addition 28 .

Conclusion
In this systematic evolutionary search for the high density energy material N 5 ASF 6 , we were able to unravel its elusive crystal structure and the cause of its marginal stability. Our investigation confirms the polymeric nitrogen configuration in this crystal to be of C 2V symmetry. The center of the Brillouin zone exhibits no negative vibrational modes and is in good agreement with the spectroscopic data. A partial softening of two of the low frequency acoustic phonon modes takes place at the vicinity of the center from different directions inside the Brillouin zone. The collective displacement of the two molecular units of this crystal in the same direction are responsible for energy transfer to higher frequency phonon, thus the explosive instability of the crystal. NBO theory also unraveled the true resonant N + 5 structure that explains its stability through an electronic encapsulation by strong triple bonding on both ends of this polymeric nitrogen. Over three decades of polymeric nitrogen investigations, many compounds were hypothesized to exist in low and high pressure domains. However, the traditional routes to find these compounds relied on intuition and trial and error to assess the possibility of their existence. In this work, evolutionary genetic predictions prove to be more efficient and require little human interaction. Moreover, PBEsol functionals within DFT is proven to be computationally more suitable than other functionals for providing a more accurate description of polymeric nitrogen structures. This work is thus believed to be a forerunner for more implementation of PBEsol in dealing with polymeric nitrogen  Table 3.
compounds in the future. The approach implemented is also assessed to successfully predict other high energy density materials, such as extended carbon-oxygen systems.

Methods
USPEX (version 9.4.4), a powerful evolutionary algorithm for stable structure predictions 29-31 starting from only their stocheometry was used. For N 5 AsF 6 under investigation, an initial pool size of 30 random crystal structures was created for the first generation. In subsequent generations, 50% of the structures were produced by heredity, 20% randomly and the remaining 30% were equally produced by permutation, soft mutation and lattice mutation. Geometry optimization for each structure was conducted in four steps, the first two steps allowing only the ions to move without changing the unit cell. In the last 2 steps, ions and the containing cell were both allowed to vary during the relaxation. Quantum ESPRESSO 32,33 , was used as the ab initio density functional theory (DFT 34,35 ) engine to carry all the structural optimizations, the electronic structure and the phonon calculations. Different DFT functionals were tried for the evolutionary search including LDA 35 and GGA 36 . Only the Perdew-Burke-Ernzerhof (PBE) functional in the framework of PBEsol 23 GGA satisfied the search criteria. Other functionals often led to the halt of the evolutional methodology of USPEX because the majority of crystal structures did not handle the system constraints. The strictest resolution in the fourth relaxation step during USPEX search had a resolution of 2π× 0.10 −1 and a kinetic energy cutoff of 680 eV. A higher cutoff of 1088 eV was imposed after the best structure was obtained in the remaining calculations. Further geometry optimization was carried out using LDA, PBE and PBEsol to assess the impact of the approximation used on the compound lattice constants. The vibrational frequencies by computing the IR activity was carried out at the PBE level of theory but with a norm conserving pseudopotential at Γ. The choice is due to the limitation in Quantum espresso's implementation of of the DFT perturbation theory in evaluating the normal modes along with their  Table 3. Comparative illustration of bond lengths and bond angles of the N + 5 between the N 5 AsF 6 found in this investigation and that of the only N + 5 salt known from the work of Vij et al 12 . The natural charge as obtained from Natural Population Analysis (NBO) is compared to the charge for the N + 5 in the work of Bartlett et al. 28 in the ideal gas phase C 2v symmetry at the the NBO(B3LYP/aug-cc-pVDZ) level of theory 28 .
Raman and IR intensities. However, the full phonon calculations were conducted using density functional perturbation theory 37 with PBEsol. The electronic Brillouin zone integration as well as phonon calculations were performed on a Monkhorst-Pack 12×12×12 k-point meshes and a 2×2×2 mesh of phonon vectors for the dynamical matrix computation. The coordinates of the two molecular units making the crystal were fed to GamessUS 38 linked with nbo.6 39 to carry out the natural bonding orbital (NBO) theory analysis without boundary conditions and without any further geometry optimization. Partial charges were then extracted from this analysis. The true bonding within N + 5 was also extracted from the bond order obtained with this NBO analysis. The NBO investigation was performed at the PBE level of theory with PBEsol and B3LYP/6-31++G. Xcrysden 40 , VESTA 41 and Jmol 42 were used for the visualization of the results and the plots of this manuscript.