A Comparative Density Functional Theory and Density Functional Tight Binding Study of Phases of Nitrogen Including a High Energy Density Material N 8

We present a comparative dispersion-corrected Density Functional Theory (DFT) and Density Functional Tight Binding (DFTB-D) study of several phases of nitrogen, including the well-known alpha, beta, and gamma phases as well as recently discovered highly energetic phases: covalently bound cubic gauche (cg) nitrogen and molecular (vdW-bound) N8 crystals. Among several tested parametrizations of N–N interactions for DFTB, we identify only one that is suitable for modeling of all these phases. This work therefore establishes the applicability of DFTB-D to studies of phases, including highly metastable phases, of nitrogen, which will be of great use for modelling of dynamics of reactions involving these phases, which may not be practical with DFT due to large required space and time scales. We also derive a dispersion-corrected DFT (DFT-D) setup (atom-centered basis parameters and Grimme dispersion parameters) tuned for accurate description simultaneously of several nitrogen allotropes including covalently and vdW-bound crystals and including high-energy phases.


Introduction
The development of new high energy density materials (HEDM) is important for a range of applications including explosives and chemical energy storage [1].Desired properties for such materials include high energy density and stability in common environments as well as under stimuli.Also desired are low cost of inputs and of synthesis as well as low environmental impact of both inputs and the products of decomposition.New HEDMs are actively researched experimentally as well as computationally [2][3][4][5][6][7][8][9].Most effective molecular HEDM contain nitrogen-and oxygen-containing heterocycles and/or NOx moieties [10].The decomposition of such molecules usually leads to the release of COx, NOx and cyanogen molecules which are toxic or environmentally damaging [2][3][4][5][6][7][8][9].
HEDMs based solely on nitrogen are attractive due to the potential for both high energy density and environmental friendliness.Nitrogen gas, N2, that is part of the atmosphere and is much used in research laboratories, is inert because the N≡N triple bond is one of the most stable chemical bonds.Alpha, beta, and gamma crystalline phases of nitrogen (molecular crystals) are well-studied.Nitrogen exhibits a uniquely large difference in energy between a single bond and a triple bond as compared to other common elements such as carbon.A large amount of energy, therefore, is released under the transformation from the single-bonded nitrogen to the molecular state.To put this statement into perspective, the average bond energy of an N-N single bond (38 kcal/mol) is considerably less than of that of a double bond (100 kcal/mol), which is around 2/5 of that of a triple bond (226 kcal/mol) [11].In comparison, the average bond energy of a C-C single bond is much higher (80 kcal/mol) compared to an N-N single bond (38 kcal/mol), while the average bond energy of C≡C (200 kcal/mol) is considerably smaller than that of a N≡N (226 kcal/mol) [11].This is to say, on decomposing from a single bonded to a triple bonded form, nitrogen would release much more energy as compared to other elements.
Singly-bonded allotropes of nitrogen are known.For example, cg-N (cubic gauche nitrogen) is a covalently bound crystalline material that had been computationally predicted [12][13][14][15] before being synthesized [16].cg-N possesses an extremely high energy density of 6.7 kcal/g [16], more than six times that of trinitrotoluene, TNT (1.1 kcal/g).However, it was not stable at room temperature and pressures below 25 GPa [16].The crystal structure of cg-N is shown in Figure 1.The discovery of cg-N is a proof simultaneously of the power of ab initio materials design (usually at the DFT, or Density Functional Theory, level [17]) and of the possibility of synthesis of nitrogen-only HEDM.However, no singly bonded allotrope of nitrogen stable at normal conditions have yet been made.
Polynitrogen compounds consisting of Nn molecular units involving single and double N-N bonds, and corresponding molecular solids stable at ambient temperatures and pressure, are expected to be easier achievable than covalent crystals.Indeed, many molecules (including all common fuels) are synthesized in a metastable state (positive heat of formation), and the vdW forces by which molecules are held together in a molecular solid are always attractive at long range.However, intermolecular interactions between Nn moieties could potentially result in decomposition [18].As of today, no molecular Nn solid has been produced, although isolated Nn molecules have been detected with n = 3-5 [19][20][21][22][23]. Recently, Hirshberg, Gerber and Krylov [24] predicted, by ab initio simulations, the existence of N8, as a (meta)stable molecular crystal, at ambient pressures.The crystal structure is shown in Figure 2. It is a vdW-bound crystal consisting of two kinds of N8 isomers, dubbed EEE and EZE in [24].The NBO (natural bond orbital) analysis of N8 allowed identifying multiple single bonds, which are responsible for the high degree of metastability [24]:  As a result, this material would release 260 kcal/mol (or 32.5 kcal per mol of N) when decomposing into N2 gas [24].This is about 4.6 kcal/g, and more than four times the gravimetric energy density of TNT and about a third smaller than that of cg-N.The lowest computed barrier to decomposition of the N8 isomers was 22 kcal/mol (0.95 eV) and the cohesive energy of the crystal 9.8 kcal/mol (0.42 eV).Clearly, this is an extremely promising HEDM.
Importantly, these metastable phases of nitrogen have been discovered computationally using DFT, and while cg-N has since been synthesized, N8 has not.DFT, while well suited to the study of bulk materials (as unit cells can be considered, which typically involve small length scales and small numbers of atoms), suffers from significant computational cost, which makes it impractical to study interfaces of such materials or the dynamics of their decomposition or reactions (e.g., of energy release).This would involve length and time scales which may be too computationally expensive.Indeed, the scaling of DFT is near-cubic with the system size (number of explicitly treated electrons or pairs thereof), which limits routinely doable calculations to simulation cells containing 10 1 -10 2 atoms.
Density Functional Tight Binding (DFTB) is an approximate DFT-based method that is about three orders of magnitude faster compared to DFT.It can be considered as a parametrization/tabulation of the most expensive parts of DFT [26][27][28][29].It achieves DFT-like accuracy for systems for which it is well parametrized.Several parameter sets have been benchmarked for DFTB with which DFTB provides ab initio accuracy for several classes of materials, including organic molecules, selected inorganic solids, and interfaces [29][30][31][32][33][34].While several parameterizations of N-N interactions for DFTB are available, the suitability of these parametrizations for the modeling of nitrogen allotropes has, however, not been established.The main purpose of this work is therefore to establish the applicability and accuracy of DFTB and its different parameterizations to the modeling of different phases of nitrogen covalent and molecular crystals including highly metastable phases like cg-N or the N8.Dispersion-corrected DFTB (DFTB-D) calculations are compared against dispersion-corrected DFT (DFT-D) calculations using atom centered basis sets.To ensure reliability of DFT-D calculations, the basis and dispersion parameters are tuned here to reproduce the structures and energetics of different phases of nitrogen.Another purpose of this work is, therefore, a DFT-D setup based on localized basis functions that is tuned for accurate description simultaneously of several nitrogen allotropes including covalently and vdW-bound crystals and including high-energy phases.The use of localized basis functions is advantageous for speed of calculations as well as for parallelizability and scalability, as they allow implementing order-N approaches.
The paper is organized as follows: Section 2 details the methodologies used in this work, Section 3 presents the results of simulations, and Section 4 concludes.

DFT Simulations
DFT simulations were performed with the SIESTA code [35] using the PBE exchange-correlation functional [36].Core electrons were simulated with Troullier-Martins pseudopotentials [37] provided with SIESTA (the 2007 version of the N.psf pseudopotentials file was used, see Appendix for pseudopotential parameters).A cutoff of 200 Ry was used for the Fourier expansion of the density, and Brillouin-zone integrations were done with a k-point Monkhorst-Pack mesh [38] with a density corresponding to at least 1 k point per 30 −1 Å −1 .Smearing equivalent to an electronic temperature of 100 K was used to improve convergence.Geometries were optimized until forces on all atoms were below 0.02 eV/Å and stresses below 0.02 GPa for N8 and 0.01 GPa for all other structures.A customized DZP (double-ξ polarized) basis set was used.The dispersion interactions in molecular crystals were modeled with the Grimme dispersion correction [39] with customized parameters (see below).All calculations were done at 0 K.

DFTB Simulations
DFTB simulations were performed using the self-consistent charge density functional tight binding scheme (SCC-DFTB) [26] and the DFTB+ code [40].The following parameter sets (Slater-Koster files) were tested: 3ob-2-1 [30,31], matsci-0-3 [32,33], and pbc-0-3 [34].These were benchmarked previously to reproduce properties of various nitrogen-containing compounds but not of nitrogen-only materials.The SCC convergence criterion was set to 1 × 10 −5 .The Brillouin zone was sampled with a k-point [38] density of no less than one point per 30 −1 Å −1 .Smearing equivalent to an electronic temperature of 300 K was used to ensure convergence.We note that, although the systems studied here are homoatomic, the self-consistent charge scheme is used because in N8, the N atoms are not equivalent and carry net charges [24].Atomic positions and lattice vectors were optimized until atomic forces are below 0.01 eV/Å.Dispersion interactions were modeled with the UFF (Universal Force Field) scheme [41,42].Note that different types of dispersion correction are available in SIESTA and DFTB+, but this is not critical as long as we achieve our purpose of reproducing structures and relative energies of different phases with both methods.

Tuning of DFT Parameters
We tuned DFT computational parameters to reproduce structural and energetic data of several nitrogen-based materials, including the N2 molecule and covalently-and vdW-bound crystals.First, the basis was tuned by adjusting the PAO.EnergyShift option in SIESTA (which effectively tunes the width of all radial ξ components of the basis functions) to reproduce the properties of the covalently bound N2 molecule and the cg-N crystal (Figure 1).The basis results in the computed bond length of N2 of 1.12 Å, in good agreement with the experimental value of 1.10 Å [43], and in the lattice constant of cg-N of 3.82 Å, in good agreement with a previously reported (zero pressure) values of 3.77-3.79Å computed with GGA (generalized gradient approximation) DFT [44,45].The experimental zero pressure estimate is 3.75 Å [16].The cg-N is 1.41 eV/atom less stable than alpha nitrogen (Table 1), in good agreement with the value of about 1.41 eV/atom reported in [24]; [46] reported a value around 1.6 eV/atom (as α-N2 is necessarily computed using a dispersion correction, this value is obtained by comparing both the energy of α-N2 obtained with the dispersion correction tuned below and the energy of cg-N computed without the dispersion correction with the energy of a free N2 molecule which is practically unaffected by the correction).The resulting PAO.EnergyShift = 0.001 Ry, and the basis set details are given in the Appendix.After the inclusion of Grimme dispersion correction with the parameters tuned as described below, the lattice constant of cg-N becomes 3.80 Å and the destabilization vs. α-N2, 1.24 eV/atom.That adding the Grimme correction makes the relative energy of cg-N worse is not surprising given the facts that by design we calibrated the basis on non-vdW systems and that due to different interatomic separations in different phases, values from different parts of the Grimme potential curve are used [47].
Table 1.Energies of different phases (eV/atom) of nitrogen vs. α -N2 obtained with DFT(-D) in SIESTA and with DFTB-D in DFTB+ with different Slater-Koster parameters.Reference plane wave DFT-D results from [48] are also given for α, β, and γ phases.The reference values for cg-N and N8 are from [24].For the α phase, absolute energies are given.Note that absolute energies are not comparable between DFT-D and DFTB-D due to differences in methods and pseudopotentials.After the basis set was tuned as described above, the Grimme parameters s6 (the scaling factor) and rij (the sum of vdW radii) were tuned to reproduce the structures and relative energies of molecular crystals α and γ nitrogen (Figure 1).The value of the C6 parameter was kept at 12.77 eV/Å 6 .These two phases are stable (local minima) at 0 pressure with the transition pressure from α (cubic phase) to γ (tetragonal phase) nitrogen of about 0.35 GPa [49].We obtained a transition pressure of 0.45 GPa.The resulting values were s6 = 2.70, rij = 4.0 Å.
These basis and Grimme parameters resulted in lattice parameters of a = 5.49 Å for α nitrogen, a = 3.85, c = 7.03 Å for β, and a = 4.03, c = 5.16 Å for γ nitrogen, in good agreement with available data [50].For example, [24] reported for α-N2 a = 5.49 Å computed with a DFT-D approach vs. the experimental value (zero pressure) of 5.65 Å [50].The lattice parameters reported in [48] (obtained with plane wave DFT calculations) are a = 5.86-5.87Å for α nitrogen, a = 4.04, c = 6.67 Å for β nitrogen, and a = 4.06, c = 5.34 Å for γ nitrogen.The experimental values are a = 4.05, c = 6.60 Å for β nitrogen (zero pressure) and a = 3.96, c = 5.11 Å for γ nitrogen (reported at 0.4 GPa) [50].This level of agreement is typical for (weakly) vdW-bound systems.The resulting energies per atom of α, β and γ nitrogen are listed in Table 1 and the lattice parameters in Table 2.The differences in the energy per atom of the phases agree very well (within typical DFT accuracy) with reference plane wave calculations (see Table 1) and show the same phase ordering (note that total energies are not comparable due to different pseudopotential approximations).
The resulting tuned DFT setup was then used to optimize the N8 crystal, which is shown in Figure 2. The lattice parameters of a = 10.66,b = 4.53, c = 6.83Å, and angles (not to be confused with the names of the phases) α = 86.5, β = 80.5, γ = 48.2degrees are in good agreement with those of [24] which reported lattice parameters within ranges a = 10.59-10.83,b = 4.40-4.67,c = 6.42-7.05Å, and α = 88.4-91.5,  = 81.7-84.9, = 44.1-47.4,depending on the DFT-D scheme used.The N8 phase is computed with the present setup to be by about 0.84 eV/atom less stable than α-N2, in good agreement with the value of about 0.88 eV/atom following from [24].
We therefore conclude that the setup derived here based on atom-centered basis functions and customized Grimme dispersion parameters (i) results in a similar level of accuracy as calculations which were reported previously for selected phases of nitrogen and (ii) provides good accuracy simultaneously for several phases, including covalent and vdW crystals and including highly metastable phases.

Comparison of DFT-D and DFTB-D Models of Different Phases of Nitrogen
After identifying basis and Grimme parameters that result in reasonably good structures and energies of several allotropes of nitrogen, we used them to benchmark several DFTB parameterizations, namely, the 3ob-2-1 [30,31], matsci-0-3 [32,33], and pbc-0-3 [34].Table 1 also shows the energy per atom of the four phases resulting from DFTB-D calculations with different parameter sets.The differences among the phases of molecular crystals agree with reference values and with our DFT-D values within typical ab initio accuracy.The lattice parameters were: a = 5.54/5.54/5.52Å for α nitrogen, a = 3.59/3.59/3.55Å, c =7.65/7.70/7.70Å for β nitrogen, and a = 4.04/4.04/4.03Å, c = 5.17/5.17/5.15Å for γ nitrogen with 3ob-2-1/matsci-0-3/pbc-0-3 DFTB parametrizations, respectively.These are summarized in Table 2 where they are compared for all phases and methods.All three parameter sets therefore appear to be suitable to model the structures of α and γ nitrogen, with a similar accuracy to that of DFT-D, and reproduce well the relative energies of α, β and γ phases.All three parameter sets have significant errors in the structure of the β phase.
For the covalently bound high-energy cg-N phase, there are significant differences in performance among the DFTB parameterizations.The 3ob-2-1 set did not result in a stable cg-N phase.When enabling the DFTB3 capability, both cg-N and N8 structures fell apart.The matsi-0-3 parameter set did result in a stable structure with a lattice constant of 3.61 Å, but its energy per atom was lower than in α-N2 (listed in Table 1).The pbc-0-3 parameter set resulted in a cg-N lattice constant of 3.73 Å, i.e., matching well the experimental estimate of 3.75 Å [16], and in destabilization vs. alpha nitrogen by 1.68 eV, which is within the ballpark of reported values (e.g., [46] reported a value around 1.6 eV; [24] reported a value of 1.41 eV/atom).
Significant differences among the three parameter sets are also observed when simulating the N8 crystal with DFTB-D.The resulting crystal structures are shown in Figure 2. Differences can be appreciated visually.Especially with 3ob-2-1, the differences in structure are drastic: there is distortion and disintegration of some N8 units.As 3ob-2-1 also did not reproduce the cg-N structure, it can be disqualified even though no significant differences in structure or energetics were noted for the low-energy phases.Both matsci-0-3 and pbc-0-3 sets reproduce the structure with similar accuracy (Table 2).Because the crystal is vdW-(i.e., weakly) bound, which can easily result in large changes in mutual orientation of molecular units with small changes of the interaction potential, the differences between the DFTB-D structures (obtained with matsci-0-3 and pbc-0-3) and the DFT-D structure (which also used a different way to account for dispersion interactions), as well as differences between matsci-0-3 and pbc-0-3 are not necessarily indicative of suitability or non-suitability of each computational setup.
Table 1 lists energy per atom differences between α nitrogen and N8 following from the three DFTB-D calculations and compares them to the DFT-D results.It is clear that the pbc-0-3 parametrization results in good agreement with DFT-D while matsci-0-3 does not.We therefore conclude that the pbc-0-3 parametrization is the only suitable to study crystal phases of nitrogen both covalently and vdW-bound and both low energy and highly metastable phases.The N8 lattice parameters with pbc-0-3 are a = 10.40,b = 4.12, c =7.76 Å, with angles α = 101.9°,β = 88.5°,γ = 49.7°.

Conclusions
We have performed dispersion-corrected Density Functional Theory and Density Functional Tight Binding simulations of different phases of nitrogen including the recently predicted N8 allotrope, which is a promising high energy density material with predicted energy release more than four times that of TNT.We tuned a DFT-D setup with atom-centered basis functions and Grimme dispersion correction to properly model covalently-as well as vdW-bound phases, including low-energy α, β, γ, nitrogen and highly metastable cg-nitrogen and N8.
DFTB is a very enticing alternative to model organic materials and specifically HEDMs, as it offers a three orders of magnitude CPU cost advantage over DFT.It is therefore important to understand if existing parameterizations of DFTB are suitable for modeling of allotropes of nitrogen, including high-energy phases.We have compared the performance of different DFTB(-D) parametrizations (3ob-2-1, matsci-0-3, and pbc-0-3) when modeling covalent and molecular crystals of different nitrogen allotropes.We conclude that only the pbc-0-3 parameter set gives result in agreement with DFT for both low-energy (α, β, γ, nitrogen) and high-energy (cg-N and N8) allotropes.This is the first time that the performance of DFTB(-D) has been benchmarked for nitrogen allotropes.

Figure 2 .
Figure 2. The unit cell of N8 as computed by (top to bottom): DFT-D and DFTB-D with the parameter sets 3ob-2-1, matsci-0-3, and pbc-0-3.Projections along (left to right): a*, b*, c* axes.(N.B.As visualization completes molecular units beyond the unit cell, their number might appear to be different in some cases).

Table 2 .
Lattice parameters of different phases nitrogen obtained with DFT-D in SIESTA and with DFTB-D in DFTB+ with different Slater-Koster parameters.The lengths are in Å and angles in degrees.