All-nitrogen anions have attracted extensive attention because of their unique chemical and physical properties, and potential applications as high-energy density materials. Here, we discovered two N44- anions with planar zigzag and tripodal structures in BeN2 compound using first-principle calculations with structural search. At ambient conditions, both of them have high kinetic and thermodynamic stability with high energy density (2.48 kJ/g, and 2.95 kJ/g relative to Be3N2 and N2 gas). The zigzag N44− anion has P21/c space group in BeN2 compound and is energetically favorite in the pressure range of 0−48 GPa. While tripodal N44− anion has Cm space group and is energetically favorite in the pressure range of 48−100 GPa. Further analysis of chemical bonding pattern and electronic properties show that the Be atoms provide 2s electrons to alter the bonding state of N44− as well as use empty outer shell 2p orbital to accommodate lone pair electrons of N atoms, forming coordinate bonds to stabilize the compounds. More important, the calculated formation of enthalpies indicate that the zigzag and tripodal N44− can be synthesized via compressing Be3N2 and N2 at a modest pressure. Our results provide a scheme to synthesis N44- anions and stabilize all-nitrogen anions by introducing beryllium atoms.

All-nitrogen salts are considered as suitable candidates for explosives or propellants.1 The high-energy content of this kind of compounds stems from the single and double bonds between nitrogen atoms. The average bond energies of N-N and N=N are 193 and 418 kJ/mol, respectively, are much less than the strong N≡N triple bond with a bonding energy of 954 kJ/mol. Therefore, singly or doubly bonded nitrogen clusters or ions decompose to triply bonded dinitrogen (N2) with an extraordinarily large energy release. Notably, all-nitrogen ion compounds are harmless to environment because the final product of the transformation is nitrogen gas. However, it is still a great challenge to synthesis a stable all-nitrogen ion compound because of the strong triple N≡N bond and relatively low activation energy toward decomposition.2 After the dinitrogen, N2, which was independently isolated in pure form from air in 1772 by Rutherford and Scheele and Cavendish,3 Curtius discovered the azide anion N3 until 1890,4 and relevant studies have been stalled since then. Thus, for the whole 20th century, only diazonium N22− and azide anion N3 were known to be stable all-nitrogen anions. Researchers have made a large number of theoretical predictions for various nitrogen derivatives from N3 to N13.5,6 However, the successful preparation of relevant compounds are few and far between. The molecule phenyl pentazole, with a five-membered all-nitrogen ring, which was characterized in the 1956.7,8 In 1999, Christe and co-workers reported the surprising synthesis of the pentanitrogen cation N5+ which was trapped as a salt with AsF6-. The purpose of the study was to produce a new type of rocket fuel that replaced the toxic steroid rocket fuel.9 Subsequently, the same group reported that the N5+SbF6 and N5+Sb2F11 salts are stable up to 70 °C. Thus, the N5+ cation became the third member of the stable all-nitrogen ions family.10 Recently, Hu’s and Lu’s groups successfully synthesized the cyclo-N5,11,12 respectively, which is a breakthrough achievement in this field. The research on all-nitrogen salts will directly promote the progress of high-energy-density materials (HEDMs). Moreover, the successful synthesis of related materials is expected to produce amazing developments in the fields of explosives, and propellants. Unfortunately, only few all-nitrogen salts is stable at ambient conditions and most of them is stabilized only at high pressure or in the presence of other elements that provide stabilization.13–16 For example, the pentazolate anion cyclo-N5 in (N5)6(H3O)3(NH4)4Cl is stabilized by chloride, ammonium, and hydronium. While the addition of non-energetic ions or groups results in the decrease of the mass ratio of cyclo-N5 in (N5)6(H3O)3(NH4)4Cl and furthermore affect the performance of explosives,11,12 it becomes essential to look for new all-nitrogen ions or method to stabilize them.

The neutral N4 molecule has been successfully isolated in the gas phase a decade ago.17 The N44- anion was first found in CsN with an open-chain structure at pressure higher than 44 GPa, all three N-N bonds in the N44- is 1.10, 1.10, and 1.26 Å.18 Compared with the neutral counterpart, the charged N44- anion, as predicted in the CsN crystal, is substantially stabilized by strong cation-anion interactions, which is expected to have an extended lifetime.18 Then, S. Zhu found N44- anions in the P21/c and Pbam structured CaN2, which is stable above 18 GPa. And it is also open-chain structure with N-N bonds 1.29, 1.37, and 1.29 Å.19 These N44- anions in CsN and CaN2 could not be preserved under ambient conditions. Therefore, it is urgent to design a new approach to synthesize the N44- anions that can be recovered at ambient conditions.

In this work, we designed a BeN2 salt by investigating the phase stabilities of BeN2 under various pressure ranging from 0 to 100 GPa. We utilized a combination of an unbiased structural search based on a particle-swarm optimization (PSO) algorithm20,21 and a first-principles density functional total energy calculation to explore the entire configuration space. The effectiveness of our method has been demonstrated by recent successes in predicting structures of various systems, ranging from elements to binary and ternary compounds.22–24 We choose beryllium atoms to stabilize the possible N44-anion after careful selection for four reasons: (1) Alkali earth metal can offer two out-shell s electrons to change the hybrid state between nitrogen atoms. (2) Beryllium element has the smallest atomic radius in the alkaline earth metal, thus resulting in the highest mass density of nitrogen and energy density. (3) The Be2+ cation and N44- anion can form a strong ionic bond to reduce the total energy of the system. (4) What’s more, the beryllium element can form rich delocalized chemical bonds, and we expect Be-N covalent bonds can bind with an adjacent nitrogen atom. This idea has been successfully used to design planar 2D Be-N sheets.25,26 We also analysed the bonding mode and electronic properties to prove the origin of high stability.

Structure searches for BeN2 were carried out at multiple pressures, e.g. 0, 5, 10, 20, 50, and 100 GPa at 0 K in the CALYPSO (Crystal structure AnaLYsis by Particle Swarm Optimization) methodology as implemented in CALYPSO code.20,21 Geometrical optimizations, total energy calculations and electronic structure calculations were performed in the framework of density functional theory within the generalized gradient approximation Perdew-Burke-Ernzerhof (GGA-PBE)27 in the Vienna Ab Initio Simulation Package (VASP) code.28 The projector-augmented-wave29 method was employed with the beryllium and nitrogen potentials which have 1s22s2 and 2s22p3 as valence states respectively, adopted from the VASP potential library. A cut off energy of 700 eV for the plane-wave basis set and appropriate k-point grid30 with a spacing of 0.03 Å−1 was used to ensure the convergence for total energies less than 1 meV/atom. To determine the dynamical stability, the phonon calculations were carried out by using a supercell approach with the finite displacement method31 through the Phonopy code.32 First-principles molecular dynamics (FPMD) simulations were performed at 300 K and 500 K to assess the thermal stability of the BeN2 crystal in Nosé-Hoover method.33,34 The chemical bonding mode of predicted BeN2 crystal is calculated by the Solid State Adaptive Natural Density Partitioning (SSAdNDP35), which allows the interpretation of chemical bonding in systems with translational symmetry in terms of classical lone pairs and two-center bonds, as well as multi-center delocalized bonding.

To seek out stable structure of N44−, we calculated the change of enthalpy (∆H) of BeN2 as shown in Fig. 1 (The whole enthalpy (∆H) of BeN2 in Fig. S1). The enthalpies of formation ∆H for BeN2 structures were calculated relative to the P21/c phase obtained in the structure search. We found that the P21/c structure transforms to a Cm structure above 48 GPa and then to a P63/mmc structure at 75 GPa. To simplify the presentation, hereafter we will denote the predicted structure as P21/c-BeN2, Cm-BeN2, and P63/mmc-BeN2.

FIG. 1.

Calculated enthalpies of formation of BeN2 under high pressure. Arrows point to the phase transition pressures. The enthalpies of formation is relative to P21/c structure with zigzag N44−, which has lowest energy at ambient conditions. Under the pressure range of 0-100 GPa, three phases of BeN2 and two phase transitions (BeN2-I (space group: P21/c) → BeN2-II (space group: Cm) → BeN2-III (space group: P63/mmc)) are predicted at 48 and 75 GPa.

FIG. 1.

Calculated enthalpies of formation of BeN2 under high pressure. Arrows point to the phase transition pressures. The enthalpies of formation is relative to P21/c structure with zigzag N44−, which has lowest energy at ambient conditions. Under the pressure range of 0-100 GPa, three phases of BeN2 and two phase transitions (BeN2-I (space group: P21/c) → BeN2-II (space group: Cm) → BeN2-III (space group: P63/mmc)) are predicted at 48 and 75 GPa.

Close modal

The predicted crystal structures are presented in Fig. 2, and corresponding lattice parameters are given in Table I. The P21/c-BeN2 and Cm-BeN2 are formed by Be2+ cations and N44− anions stabilizes into a monoclinic structure in Fig. 2a-b. Each two Be atoms are coordinated with four N atoms in these structures. Interestingly, we find polymerization of nitrogen is achieved by forming zigzag (Fig. 2a) or tripod nitrogen anion (Fig. 2b). Calculated N-N distances in the zigzag N44− are 1.367 Å and 1.298 Å, and these in tripod N44− are 1.263 Å and 1.381 Å, which are between the standard bond lengths for the single and double bonds (1.25 Å and 1.45 Å).36,37 In these structures, nitrogen is fully reduced to the N44− anion with complete subshells. The charged N44− as predicted in the BeN2 crystal is substantially more stable than the neutral counterpart through strong cation-anion interactions, which is expected to have an extended lifetime. Furthermore, compared with the neutral counterpart, due to the additional electrons obtained from Be, the N44− anion has a lower extent of electrons-sharing results in a higher energy content. The energy densities of P21/c-BeN2 and Cm-BeN2 are 2.48 kJ/g, and 2.95 kJ/g, respectively, relative to Be3N2 and N2 gas, which are close to cyclo-N5 in LiN5 (2.72 kJ/g).5 The predicted P63/mmc-BeN2 (Fig. 2c) is a typical pernitride compounds38–41 with Be2+ cations and N22− anions. Each Be atom is coordinated with two N atoms in this structure and the distance of N-N is 1.243 Å.

FIG. 2.

Crystalline structures of the predicted stable BeN2 phases. (a) P21/c structure of BeN2 at ambient pressure (b) Cm structure of BeN2 at ambient pressure (c) P63/mmc structure of BeN2 at ambient pressure. The smaller spheres are nitrogen atoms, and the larger ones are beryllium atoms.

FIG. 2.

Crystalline structures of the predicted stable BeN2 phases. (a) P21/c structure of BeN2 at ambient pressure (b) Cm structure of BeN2 at ambient pressure (c) P63/mmc structure of BeN2 at ambient pressure. The smaller spheres are nitrogen atoms, and the larger ones are beryllium atoms.

Close modal
TABLE I.

The unit-cell parameters and atomic positions of the P21/c, Cm, and P63/mmc phase at 0GPa, respectively.

Space groupLattice parameter (Å)AtomWyckoff positionxyz
P21/c a= 3.8518 α= 90.0000 Be 4e 0.901130 0.341930 0.378940 
 b= 4.4387 β= 70.5265 N1 4e 0.155490 0.222460 0.087330 
 c= 5.7054 γ= 90.0000 N2 4e 0.493840 0.354540 0.990490 
Cm a= 4.5503 α= 90.0000 Be 4b 0.825880 0.830710 0.856110 
 b= 7.8929 β= 121.1099 N3 4b 0.557280 0.850750 0.103210 
 c= 2.8978 γ= 90.0000 N1 2a 0.086660 0.000000 0.980620 
   N2 2a 0.382350 0.000000 0.041650 
P63/mmc a=2.6631 α=90.0000 Be 2a 0.000000 0.000000 0.000000 
 b=2.6631 β=90.0000 4f 0.666670 0.333333 0.159490 
 c=6.8677 γ=120.0000      
Space groupLattice parameter (Å)AtomWyckoff positionxyz
P21/c a= 3.8518 α= 90.0000 Be 4e 0.901130 0.341930 0.378940 
 b= 4.4387 β= 70.5265 N1 4e 0.155490 0.222460 0.087330 
 c= 5.7054 γ= 90.0000 N2 4e 0.493840 0.354540 0.990490 
Cm a= 4.5503 α= 90.0000 Be 4b 0.825880 0.830710 0.856110 
 b= 7.8929 β= 121.1099 N3 4b 0.557280 0.850750 0.103210 
 c= 2.8978 γ= 90.0000 N1 2a 0.086660 0.000000 0.980620 
   N2 2a 0.382350 0.000000 0.041650 
P63/mmc a=2.6631 α=90.0000 Be 2a 0.000000 0.000000 0.000000 
 b=2.6631 β=90.0000 4f 0.666670 0.333333 0.159490 
 c=6.8677 γ=120.0000      

Structures determined by geometric optimization may be affected by kinetic instabilities. To examine the stability of the P21/c, Cm, and P63/mmc structures, we calculate the phonon spectra. Their phonon spectra and density of states (PHDOS) of the P21/c, Cm, and P63/mmc structures at 0 GPa are shown in Fig. 3. For P21/c-BeN2 and Cm-BeN2, no imaginary vibrational modes over the Brillouin (Fig. 3a-b) confirm their dynamical stabilities. However, it is with regret that P63/mmc-BeN2 has seriously imaginary vibrational modes and even at high pressure (Fig. 3c, S8, and S9), which indicates that pernitride in BeN2 is unstable and different from the heavy alkaline-earth metals pernitride.38–44 Hence, we will focus on the P21/c-BeN2 and Cm-BeN2 in this work. On the basis of the analysis of the phonon density of states of P21/c (Fig. 3a) and Cm structure (Fig. 3b), the low frequency vibration mode mainly comes from the strong coupled vibration between Be and N atoms, while the high frequency vibration is mainly attributed to NN stretching mode.

FIG. 3.

Phonon dispersions and PHDOS of BeN2. Phonon dispersions and PHDOS of (a) P21/c-BeN2, (b) Cm-BeN2, and (c) P63/mmc-BeN2 at 0 GPa, where the blue and red lines indicate the PHDOS of Be and N atoms, respectively.

FIG. 3.

Phonon dispersions and PHDOS of BeN2. Phonon dispersions and PHDOS of (a) P21/c-BeN2, (b) Cm-BeN2, and (c) P63/mmc-BeN2 at 0 GPa, where the blue and red lines indicate the PHDOS of Be and N atoms, respectively.

Close modal

To further confirm the thermodynamic stability at atmospheric pressure, we performed first-principles molecular dynamics simulations (FPMD) using a 2 x 2 x 2 supercell at the temperature 300 K with a time step of 1 fs. As shown in Fig. 4a and 4c, the image of the geometry structure at the end of the 10 ps simulations clearly shows that both zigzag and tripodal N44- anion maintain their structural integrities, except for some thermal fluctuations. The potential energy of the analog system fluctuates around 1522 eV and 1508 eV (Fig. 4b and 4d), respectively. Obviously, the P21/c and Cm phases are thermodynamically stable at ambient conditions and might be recoverable.

FIG. 4.

Molecular Dynamics Simulation of BeN2 at 0 GPa. The equilibrium structures of (a) P21/c and (c) Cm at the end of 10 ps FPMD simulations at the temperatures of 300 K. (b) and (d) the fluctuations of the total potential energies of P21/c and Cm structures with respect to FPMD steps at 300K, respectively.

FIG. 4.

Molecular Dynamics Simulation of BeN2 at 0 GPa. The equilibrium structures of (a) P21/c and (c) Cm at the end of 10 ps FPMD simulations at the temperatures of 300 K. (b) and (d) the fluctuations of the total potential energies of P21/c and Cm structures with respect to FPMD steps at 300K, respectively.

Close modal

We have calculated the electron band structures and the corresponding projected density of states (PDOS) to detect the electronic properties of the P21/c-BeN2 and Cm-BeN2. The calculated band structure and PDOS (Fig. 5a, 5b, and 5c) show that the P21/c-BeN2 crystal is a semiconductor with a band gap of 0.79 eV. The Cm-BeN2 crystal is also a semiconductor with a band gap of 1.64 eV (Fig. 5e, 5f, and 5g). Since density functional calculations often result in a fairly underestimation of the energy gap, it is expected that the actual band gap will be larger than the calculated results. Furthermore, there are obvious overlap between N 2p and Be 2p for both P21/c-BeN2 (Fig. 5b and 5c) and Cm-BeN2 (Fig. 5f and 5g), indicating strong coupling between Be 2p and N 2p orbitals.

FIG. 5.

Electronic properties of the P21/c and Cm structures of BeN2 at 0 GPa. (a) Band structure and projected density of states for (b) Be and (c) N atoms. (d) Crystal orbital Hamilton populations (COHP) between Be and N atoms of P21/c; (e) Band structure and projected density of states for (f) Be and (g) N atoms. (h) COHP between Be and N atoms of Cm.

FIG. 5.

Electronic properties of the P21/c and Cm structures of BeN2 at 0 GPa. (a) Band structure and projected density of states for (b) Be and (c) N atoms. (d) Crystal orbital Hamilton populations (COHP) between Be and N atoms of P21/c; (e) Band structure and projected density of states for (f) Be and (g) N atoms. (h) COHP between Be and N atoms of Cm.

Close modal

To clarify the bonding case, the crystal orbital Hamiltonian group (COHP) between Be and N atoms is calculated in the LOBSTER program.45 COHP is a measure of the overlap strength and an energy resolved quantity of crystal orbital overlap population which is negative for bonding states and positive for antibonding states.46 The calculated COHP between Be and nearest N atoms, as shown in Fig. 5d and 5h, demonstrates the strong covalent bond feature between them in both P21/c-BeN2 and Cm-BeN2.

To further describe the bonding feature of N44- anions in P21/c and Cm phases, we use the SSAdNDP software to analyze the chemical bonding patterns. According to the chemical bonding analysis (Fig. 6), in a primitive cell, there are sixteen Be-N σ bonds (Fig. 6a and 6e) which require 2 × 16 = 32 electrons; as well as six 2c-2e N-N σ bonds (Fig. 6b and 6f) and two 2c-2e π bonds (Fig. 6c and 6g), which requires 2 × 6 + 2 × 2 = 16 electrons. The sum of the electrons is 48, This is in completely agreement with the number of outer valence electrons in the BeN2 unit cell (one unit cell contains two BeN2 formulations, 2s2 and 2s22p3 for Be and N elements, respectively). The Occupation numbers of the bonds are close to the ideal value, so the reliability of the chemical bonding picture is high. In the zigzag N44- anion of the P21/c-BeN2 crystal, the terminal atoms are bonded to the internal atoms by a single bond (σ) and the two internal atoms are double bonded by an N=N bond (σ and π). At the same time, each terminal N atom is bonded to three nearest Be atoms and each internal N atom is bonded to one nearest Be atom via Be-N σ bonds (Fig. 6d). For tripodal N44- anion in Cm-BeN2, one of the terminal N atoms is bonded to the central atom by an N=N double bond (σ and π) and the other two atoms are bonded to the central atom by a N-N single bond (σ). In the meantime, the N atom with double bond is bonded to two nearest Be atoms and the N atom with single bond is bonded to three nearest Be atoms via Be-N σ bonds (Fig. 6h).

FIG. 6.

Chemical bonding pattern. Schematic of the SSAdNDP chemical bonding pattern for P21/c - BeN2 and Cm- BeN2 unit cell with (a) and (e) 16 × 2c-2e Be-N σ bonds, occupation number threshold: 1.70 |e|, isovalue = 0.1; (b) and (f) 6 × 2c-2e N-N σ bonds, occupation number threshold: 1.70 |e|, isovalue = 0.1 and 0.2, respectively; (c) and (g) 2 × 2c-2e N-N π bonds, occupation number threshold: 1.70 |e|, isovalue = 0.1. (d) and (h) Bonding state between N and N atoms, and between N and Be atoms for P21/c - BeN2 and Cm- BeN2.

FIG. 6.

Chemical bonding pattern. Schematic of the SSAdNDP chemical bonding pattern for P21/c - BeN2 and Cm- BeN2 unit cell with (a) and (e) 16 × 2c-2e Be-N σ bonds, occupation number threshold: 1.70 |e|, isovalue = 0.1; (b) and (f) 6 × 2c-2e N-N σ bonds, occupation number threshold: 1.70 |e|, isovalue = 0.1 and 0.2, respectively; (c) and (g) 2 × 2c-2e N-N π bonds, occupation number threshold: 1.70 |e|, isovalue = 0.1. (d) and (h) Bonding state between N and N atoms, and between N and Be atoms for P21/c - BeN2 and Cm- BeN2.

Close modal

Both zigzag and tripodal N44- anions in BeN2 exhibit high kinetic and thermodynamic stabilities due to the existence of two kind interactions between Be2+ cations and N44- anions. One is Coulomb interactions. Two beryllium atoms lose the 2s electrons and became Be2+ cations in BeN2 crystal. At the same time, the N4 gets four electrons and becomes N44-. Thus, a Coulomb interaction between Be2+ and N44- is formed. The other is covalent bond. All the chemical bonding pattern analyses, PDOS and COHP calculations, illustrate the covalent bonding features between Be2+ and N44- anion. In one word, the combination of the Coulomb interactions and covalent bonds between Be2+ cations and N44- ions are collectively responsible for the high stabilities of zigzag and tripodal N44- anions in BeN2 crystal.

The beryllium nitride, Be3N2, with two known phases (alpha and beta) that has been demonstrated to be synthesizable in experiment. The alpha phase is a cubic bcc structure, which is stable between 20–1200 °C47,48 and at temperatures over 1400 °C changes to a hexagonal structure, which is the beta phase.49 To investigate a possible high-pressure experimental rout to P21/c-BeN2 and Cm-BeN2, the enthalpies of formation of two predicted BeN2 relative to Be3N2 and N2 have been calculated as shown in Fig. 7. The results show that the P21/c-BeN2 has the lowest energy in the pressure range of 22-48 GPa and the Cm-BeN2 become the energetically favorite one above 48 GPa. Which indicate the zigzag N44− can be synthesized in the pressure range 22-48 GPa and the tripodal N44− can be synthesized in the pressure range 48-75 GPa considering that the unstable P63/mmc-BeN2 has lowest energy above 75 GPa via compressing Be3N2 and N2.

FIG. 7.

Enthalpies of formation of P21/c-BeN2 and Cm-BeN2 relative to Be3N2 and N2.

FIG. 7.

Enthalpies of formation of P21/c-BeN2 and Cm-BeN2 relative to Be3N2 and N2.

Close modal

We have reported two all-nitrogen anion N44− in the BeN2 compounds based on thorough structural searches and density functional calculations. Both zigzag and tripodal N44− in the BeN2 have high kinetic and thermodynamic stabilities with high energy densities (2.48 and 2.95 kJ/g) at ambient conditions. The beryllium atoms in BeN2 serve as electron donors which lose the 2s electrons and become Be2+. Thus, a coulomb interaction between Be2+ cations and N44- anions is formed. The analyses of electronic properties and chemical bonding pattern demonstrate the covalent bonds between Be2+ and adjacent N also enhance the stability of BeN2 crystal. The calculated formation of enthalpies of P21/c-BeN2 and Cm-BeN2 indicate that the zigzag and tripodal N44− can be synthesized via compressing Be3N2 and N2 at a modest pressure. We expect that our theoretical research could encourage the experimental realization of all-nitrogen N44− anions.

See supplementary material for the complete enthalpies of BeN2 and phonon dispersion curves of the predicted structures at high pressure.

This work was supported by the National Natural Science Foundation of China under Grants Nos. 11674144 and 11774128, the Natural Science Foundation of Shandong Province Nos. JQ201602, ZR2018MA038, and ZR2018JL003.

1.
G.
Steinhauser
and
T. M.
Klapötke
,
Angew. Chemie Int. Ed.
47
,
3330
(
2008
).
2.
K. O.
Christe
,
Science
355
,
351
(
2017
).
3.
N.
Greenwood
and
A.
Earnshaw
,
Chemistry of the Elements
(
Elsevier
,
1997
).
4.
T.
Curtius
,
Berichte Der Dtsch. Chem. Gesellschaft
23
,
3023
(
1890
).
5.
F.
Peng
,
Y.
Yao
,
H.
Liu
, and
Y.
Ma
,
J. Phys. Chem. Lett.
6
,
2363
(
2015
).
6.
A.
Kulkarni
,
J. C.
Schön
,
K.
Doll
, and
M.
Jansen
,
Chem. - An Asian J.
8
,
743
(
2013
).
7.
R.
Huisgen
and
I.
Ugi
,
Angew. Chemie
68
,
705
(
1956
).
8.
R.
Huisgen
and
I.
Ugi
,
Chem. Ber.
90
,
2914
(
1957
).
9.
K. O.
Christe
,
W. W.
Wilson
,
J. A.
Sheehy
, and
J. A.
Boatz
,
Angew. Chemie Int. Ed.
38
,
2004
(
1999
).
10.
A.
Vij
,
W. W.
Wilson
,
V.
Vij
,
F. S.
Tham
,
J. A.
Sheehy
, and
K. O.
Christe
,
J. Am. Chem. Soc.
123
,
6308
(
2001
).
11.
Y.
Xu
,
Q.
Wang
,
C.
Shen
,
Q.
Lin
,
P.
Wang
, and
M.
Lu
,
Nature
549
,
78
(
2017
).
12.
C.
Zhang
,
C.
Sun
,
B.
Hu
,
C.
Yu
, and
M.
Lu
,
Science
355
,
374
(
2017
).
13.
W.
Zhang
,
K.
Wang
,
J.
Li
,
Z.
Lin
,
S.
Song
,
S.
Huang
,
Y.
Liu
,
F.
Nie
, and
Q.
Zhang
,
Angew. Chemie Int. Ed.
57
,
2592
(
2018
).
14.
C.
Sun
,
C.
Zhang
,
C.
Jiang
,
C.
Yang
,
Y.
Du
,
Y.
Zhao
,
B.
Hu
,
Z.
Zheng
, and
K. O.
Christe
,
Nat. Commun.
9
(
1
),
1269
(
2018
).
15.
C.
Zhang
,
C.
Yang
,
B.
Hu
,
C.
Yu
,
Z.
Zheng
, and
C.
Sun
,
Angew. Chemie Int. Ed.
56
,
4512
(
2017
).
16.
Y.
Xu
,
P.
Wang
,
Q.
Lin
, and
M.
Lu
,
Dalt. Trans.
46
,
14088
(
2017
).
17.
F.
Cacace
,
G.
De Petris
, and
A.
Troiani
,
Science
295
,
480
(
2002
).
18.
F.
Peng
,
Y.
Han
,
H.
Liu
, and
Y.
Yao
,
Sci. Rep.
5
,
16902
(
2015
).
19.
S.
Zhu
,
F.
Peng
,
H.
Liu
,
A.
Majumdar
,
T.
Gao
, and
Y.
Yao
,
Inorg. Chem.
55
,
7550
(
2016
).
20.
Y.
Wang
,
J.
Lv
,
L.
Zhu
, and
Y.
Ma
,
Phys. Rev. B
82
,
094116
(
2010
).
21.
Y.
Wang
,
J.
Lv
,
L.
Zhu
, and
Y.
Ma
,
Comput. Phys. Commun.
183
,
2063
(
2012
).
22.
X.
Wang
,
Y.
Wang
,
M.
Miao
,
X.
Zhong
,
J.
Lv
,
T.
Cui
,
J.
Li
,
L.
Chen
,
C. J.
Pickard
, and
Y.
Ma
,
Phys. Rev. Lett.
109
,
175502
(
2012
).
23.
J.
Li
,
L.
Sun
,
X.
Wang
,
H.
Zhu
, and
M.
Miao
,
J. Phys. Chem. C
122
,
22339
(
2018
).
24.
X.
Wang
,
J.
Li
,
J.
Botana
,
M.
Zhang
,
H.
Zhu
,
L.
Chen
,
H.
Liu
,
T.
Cui
, and
M.
Miao
,
J. Chem. Phys.
139
,
164710
(
2013
).
25.
F.
Li
,
Y.
Wang
,
H.
Wu
,
Z.
Liu
,
U.
Aeberhard
, and
Y.
Li
,
J. Mater. Chem. C
5
,
11515
(
2017
).
26.
X.
Li
,
S.
Zhang
,
C.
Zhang
, and
Q.
Wang
,
Nanoscale
10
,
18936
(
2018
).
27.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
28.
G.
Kresse
and
J.
Furthmüller
,
Phys. Rev. B - Condens. Matter Mater. Phys.
54
,
11169
(
1996
).
29.
P. E. P.
Blöchl
,
Phys. Rev. B
50
,
17953
(
1994
).
30.
J. D.
Pack
and
H. J.
Monkhorst
,
Phys. Rev. B
16
,
1748
(
1977
).
31.
K.
Parlinski
,
Z. Q.
Li
, and
Y.
Kawazoe
,
Phys. Rev. Lett.
78
,
4063
(
1997
).
32.
A.
Togo
,
F.
Oba
, and
I.
Tanaka
,
Phys. Rev. B
78
,
134106
(
2008
).
33.
S.
Nosé
,
J. Chem. Phys.
81
,
511
(
1984
).
34.
W. G.
Hoover
,
Phys. Rev. A
31
,
1695
(
1985
).
35.
T. R.
Galeev
,
B. D.
Dunnington
,
J. R.
Schmidt
, and
A. I.
Boldyrev
,
Phys. Chem. Chem. Phys.
15
,
5022
(
2013
).
36.
Y.
Morino
,
Y.
Nakamura
, and
T.
Iijima
,
J. Chem. Phys.
32
,
643
(
1960
).
37.
M.
Carlotti
,
J. W. C.
Johns
, and
A.
Trombetti
,
Can. J. Phys.
52
,
340
(
1974
).
38.
S.
Yu
,
B.
Huang
,
Q.
Zeng
,
A. R.
Oganov
,
L.
Zhang
, and
G.
Frapper
,
J. Phys. Chem. C
121
,
11037
(
2017
).
39.
H.
Wang
,
Y.
Yao
,
Y.
Si
,
Z.
Wu
, and
G.
Vaitheeswaran
,
J. Phys. Chem. C
118
,
650
(
2014
).
40.
G.
Auffermann
,
Y.
Prots
, and
R.
Kniep
,
Angew. Chemie Int. Ed.
40
,
547
(
2001
).
41.
G. V.
Vajenine
,
G.
Auffermann
,
Y.
Prots
,
W.
Schnelle
,
R. K.
Kremer
,
A.
Simon
, and
R.
Kniep
,
Inorg. Chem.
40
,
4866
(
2001
).
42.
S. B.
Schneider
,
R.
Frankovsky
, and
W.
Schnick
,
2
(
2012
).
43.
S. B.
Schneider
,
M.
Mangstl
,
G. M.
Friederichs
,
R.
Frankovsky
,
J.
Schmedt Auf Der Günne
, and
W.
Schnick
,
Chem. Mater.
25
,
4149
(
2013
).
44.
B. R.
Brooks
,
C. L.
Brooks
 III
,
J. A. D.
Mackerell
,
L.
Nilsson
,
R. J.
Petrella
,
B.
Roux
,
Y.
Won
,
G.
Archontis
,
C.
Bartels
,
S.
Boresch
,
A.
Caflisch
,
L.
Caves
,
Q.
Cui
,
A. R.
Dinner
,
M.
Feig
,
S.
Fischer
,
J.
Gao
,
M. W. I.
Hodoscek
, and
M.
Karplus
,
J. Comput. Chem.
30
,
1545
(
2009
).
45.
S.
Maintz
,
V. L.
Deringer
,
A. L.
Tchougréeff
, and
R.
Dronskowski
,
J. Comput. Chem.
37
,
1030
(
2016
).
46.
R.
Dronskowski
and
P. E.
Bloechl
,
J. Phys. Chem.
97
,
8617
(
1993
).
47.
G. T.
Furukawa
and
M. L.
Reilly
,
74A
,
617
(
1970
).
48.
T. B.
Douglas
and
W. H.
Payne
,
J. Res. Notional Bur. Stand. - A. Phys. Chem.
73A
,
471
(
1969
).
49.
D.
Hall
,
G. E.
Gurr
, and
G. A.
Jeffrey
,
Zeitschrift Für Anorg. Und Allg. Chemie
369
,
108
(
1969
).

Supplementary Material