Antifluorite-type Na₅FeO₄ as a low-cost, environment-friendly cathode with combined cationic/anionic redox activity for sodium ion batteries: a first-principles investigation

Abstract

The rapid electrification of our society and the transition towards a larger share of intermittent renewable energy sources in our electricity grids will dramatically increase the demand for cheap energy storage. Sodium ion batteries (SIBs) show a lot of promise to provide the required stationary storage at the grid level at low cost owing to the natural abundance and geographical availability of sodium. In addition, alkali-rich cathode materials exhibiting anionic redox contributions have garnered much attention over the past decade as a strategy to increase the specific capacity. In this work, we investigate for the first time the sodium-rich compound Na₅FeO₄ as a potential low-cost, environment-friendly cathode for sodium ion batteries from first principles using density functional theory. We investigate three low-energy polymorphs related to the antifluorite structure, verify their dynamical and mechanical stabilities, and show that they exhibit promising ion diffusive properties. As alkali-rich cathode materials are prone to oxygen loss during cycling, we investigate cycling stability with respect to phase transformations and oxygen loss and identify in particular one promising cycling interval that can reversibly shuttle 1.5 Na⁺ per formula unit between Na₅FeO₄ and Na_3.5FeO₄ with a gravimetric energy density exceeding 360 W h kg⁻¹. Investigations into possible redox mechanisms reveal that the charge compensation occurs simultaneously on Fe- and O-atoms in FeO₄-tetrahedra, which suggests that Na₅FeO₄, if realised experimentally as a cathode material, would join the family of combined cationic/anionic redox compounds.

---

In the past decade, sodium ion batteries (SIBs) have attracted increasing attention as an alternative to lithium ion batteries (LIBs) for applications where the energy density is not of utmost concern, such as for stationary grid storage.¹ For these applications, the levelised cost of energy (LCOE) is the most important factor, and this will in large depend on manufacturing costs and cycle lifetime. One way to lower the cost of the electrode materials is by only utilising earth abundant elements that are not prone to price volatility. SIBs are interesting in this respect due to the high abundance, geographical availability and correspondingly low cost of sodium compared to lithium. In addition, by only using iron as a transition metal redox centre, the total material cost of the cathode could indeed be very low. Furthermore, by employing SIBs for applications where the higher energy densities found in LIBs are not strictly required, regions without a natural abundance of lithium resources would be able to lessen their import dependence on critical raw materials (CRM) required for battery production, such as within the European Union where such a CRM independence is a stated goal.² Continued concerns about lithium resource availability could also allow SIBs to play a larger role in the future energy system, alleviating pressure on lithium where it is not required and further contribute to CRM independence.^3–6

Another topic that has generated a lot of interest in the academic community over the past decade is cathode materials that exhibit combined cationic and anionic redox activity.⁷ In these materials, as Li⁺/Na⁺ is removed from the compound during charging, there is a redox compensation not only from the oxidation of the transition metal, but also from a partial oxidation of an oxygen anion. This has been shown to arise from the existence of labile oxygen states close to the Fermi level.⁸ These alkali metal rich materials provide the essential environments surrounding the oxygen anions to enable the labile states, and at the same time substitutes heavier transition metal ions for the lighter alkali metal ions resulting in a higher theoretical specific capacity. Many of the researched materials have been layered structures, such as the proof-of-concept cathode Li₂RuO₃,⁹ and while most such materials have been prone to poor cycling stability due to an irreversible loss of oxygen from the lattice, advancements in our fundamental understanding of such materials suggest that the issues of structural instability at high deintercalation levels should be less pronounced in tridimensional structures.⁸

One such class of tridimensional materials are the antifluorite-type compounds, such as Li₅FeO₄, Li₆CoO₄ and Li₆MnO₄, which have all been studied as potential cathode materials for LIBs.^10,11 Li₅FeO₄ has, for example, been shown to allow reversible extraction of up to 1.2 Li⁺, giving a specific capacity of just above 200 mA h g⁻¹.¹² The crystal structures of these compounds can all be described by considering a 2 × 2 × 2 supercell of the ideal antifluorite structure Li₂O, with oxygen anions forming a cubic closest packed sublattice with lithium cations occupying all tetrahedral holes. Comparing Li₅FeO₄ with Li₂O (Li₈O₄), three Li⁺ are replaced by one Fe³⁺ and two Li⁺-vacancies (v′_Li). In Li₆CoO₄, two Li⁺ are replaced by one Co²⁺ and one v′_Li. Their sodium analogues were reported to exist in the 70s and 80s by Hoppe and co-workers^13–16 although with slightly different stoichiometries. In these cases, the sodium contents of the Co- and Mn-structures are lower, compensated by a higher oxidation state in the transition metal cations, and exhibiting other variants of the antifluorite structure. To our knowledge, there are no reports to date of these materials having been investigated for possible use as cathode materials for SIBs.

In this work, we investigate the sodium oxoferrate Na₅FeO₄ as a potential candidate as a cheap, environmentally friendly cathode for SIBs by use of density functional theory. After describing the computational details, we show results from an investigation of the phase stability based on experimentally known and theoretically predicted crystal structures of A₅MO₄ (A = alkali metals, M = transition metals) from the structure databases ICSD¹⁷ and Materials Project,¹⁸ and show how the relative stability of the most stable modifications shift with temperature. We verify that the low-energy modifications are all dynamically and mechanically stable, and verify that their diffusive properties are suitable for use as a cathode material based on migration barriers obtained from nudged elastic band-calculations (NEB). In order to predict how much Na⁺ can be expected to be reversibly extracted from the structure, we compare the phase stability of the desodiated phases with oxygen deficient structures and calculate the formation enthalpies for oxygen vacancies at different desodiation levels. Lastly, we discuss the possible redox mechanisms during desodiation from a structural and electronic point of view.

1 Computational details

Density functional theory (DFT) calculations are performed using the projected-augmented-wave (PAW)^19,20 implementation of the Vienna Ab initio Simulation Package (VASP 5.4.4).^21–24 The Perdew, Burke and Ernzerhof (PBE) functional is used for the exchange-correlation term.²⁵ The Hubbard parameter U is used following the rotationally invariant form by Dudarev and Botton.²⁶ The effective U-value for the d-orbitals is set to 5.3 eV for Fe in correspondance with calibrated values for oxide systems by the Materials Project, which were obtained by fitting to experimental binary formation enthalpies.¹⁸ The optimised structures are obtained by minimising the total energy with a convergence criterion of 10⁻³ eV using the conjugate gradient algorithm. The convergence threshold for the electronic self-consistent calculations are set to 10⁻⁵ eV. Brillouin zone (BZ) integration is performed with a Γ-centered grid with a smallest allowed spacing between k-points of 0.4 Å⁻¹. The BZ integration is done according to the Tetrahedron method with Blöchl-corrections for construction of the energy–volume curves and electronic structure calculations.²⁷ For subsequent calculations, Guassian smearing with broadening width of 0.2 eV was used. We explicitly include all 2p and 3s electrons for Na, all 2s and 2p electrons for O and all 3p, 3d and 4s electrons for Fe. The energy cutoff is set to 500 eV, and the calculations are converged to within 2.5 meV per atom with respect to a higher energy cutoff of 800 eV. The energy–volume curves are obtained by stepwise relaxation of the structure. In the first step, only the ionic positions are relaxed, in the second step both ionic positions and cell shape are relaxed, and finally both ionic positions, cell shape and cell volume are relaxed. From the relaxed structure we then construct fixed volume configurations by applying a diagonal strain in increments of 1% volume change in both the tensile and compressive directions, and allowing ionic positions to relax. The obtained energies are then fitted to the Birch–Murnaghan equation of state²⁸ to obtain equilibrium values for the total energy, volume, bulk modulus and its pressure derivative.

The phonopy package is used to construct displacements for all phonon calculations.²⁹ For these calculations, we construct supercells when needed to ensure that each lattice parameter is at least 10 Å. First, we allow the structure to relax using a stricter convergence criterion with a force convergence threshold of 10 meV Å⁻¹ and electronic convergence criterion of 10⁻⁸ eV. The forces are then calculated for all displacements using the same electronic convergence criterion where all projection operators are evaluated in reciprocal space for added accuracy, as opposed to the automatic scheme provided by VASP used in all other calculations. We employ the finite differences method with displacements lengths of 0.1 Å. The phonopy package is further used to calculate the phonon dispersion relations, phonon density of states and thermal properties.

The single crystal elastic constants are calculated using the finite strain method, and the resulting symmetrised tensor is diagonalised using the NumPy package³⁰ to obtain the eigenvalues.

Electronic density of states have been calculated by running a static, self-consistent calculation on a pre-relaxed structure with the smallest allowed spacing between k-points reduced to 0.1 Å⁻¹ for a finer grid in k-space.

To assess the electrochemical behaviour, we calculate the formation energies ΔE_form of several Na-deficient compositions, Na_5−xFeO_4−y (x = 0.0–3.0, y = 0.0–2.0), to construct the convex hull as outlined by Urban et al.³¹ from

ΔE_form(1 − X) = E(1 − X) − (1 − X)E(1) − XE(0) (1)

where the E(1) and E(0) are the total energies of the end compositions and E(1 − X) (0 ≤ X ≤ 1) is the total energy of an intermediate composition a distance X from E(0) and (1 − X) from E(1), from which E_form(0) and E_form(1) are set to 0 by construction. We used pymatgen to obtain the Na-deficient structures by generating all combinations of vacant structures and proceeding with the set of ten with the lowest unique Ewald sum for each of the deintercalation steps.³² In addition, we performed an extensive search, using matminer to search the Materials Project database³³ and manually searching the ICSD-database, for structures matching the stoichiometry of the deintercalated structures, as well as for structures with a lower oxygen stoichiometry. The latter structures were included in order to assess when irreversible oxygen release is expected, and to compare these modifications with structures of higher oxygen content, we calculate their total energies as

where E_{Na5−xFeO4−y} is the ground state energy of the oxygen deficient modification and E_O2 is the total energy of a gas phase oxygen molecule.

We further treat the most stable desodiated structures along the convex hull and calculate the equilibrium voltages vs. Na⁺/Na(V) from the following equation:

where (x₂ − x₁) is the number of deintercalated Na⁺-ions for the given voltage step and E(x_i) is the total energy for the respective phase. The chemical potential μ_Na is set to the total energy per atom of the bcc form of Na metal.

The activation energies for the diffusion of Na⁺-ions are calculated using the NEB-method.³⁴ First initial and final states of a given jump are created, before a set of five images are generated between these states from linear interpolation between the end points using VTST Tools.^34,35 When needed, we create a supercell of minimum 10 Å × 10 Å × 10 Å to ensure sufficient spacing between the vacancy and its image in the neighbouring unit cells to avoid self-interaction in order to investigate diffusion in the dilute limit. We first conduct a NEB calculation to obtain an initial guess of the minimum energy path. Then the calculation is restarted using the climbing nudged elastic band (cNEB) method to obtain a more accurate result.³⁵ The initial vacant structures are created by considering all Na–Na distances shorter than 3 Å, in certain cases extended somewhat to complete a percolating diffusion path. After removal of the appropriate Na-ions for each jump, the atomic positions are relaxed before creating the images. The ionic optimisation for the intermediate images are carried out using the FIRE algorithm³⁶ with a force convergence threshold of 40 meV Å⁻¹. The diffusion coefficients of Na⁺ (D_Na⁺) as a function of temperature are estimated from

D_Na⁺(T) = d²ν₀exp(−E_a/k_BT) (4)

where d is the hopping distance, ν₀ is the attempt frequency (assumed to be 10 THz³⁷), E_a is the activation barrier obtained by the NEB calculations and k_B is the Boltzmann constant.

The oxygen vacancy formation enthalpies are calculated based on the following reaction:

where 0 < x < 3 and δ is dependent on the size of the unit cell used, but in all cases limited to give a vacancy concentration below 3.13%. The vacancy formation energies ΔE_vac are calculated as

where E_Na5−xFeO4, E_{Na5−xFeO4−δ} and E_O2 are the total energies of the pristine structures, the oxygen deficient structures and gas phase oxygen molecule, respectively. The oxygen deficient structures are generated in a similar manner as described above for the Na-deficient structures using pymatgen. To study the thermodynamic stability of the oxygen vacancies as a function of temperature, we substitute the total energy of the oxygen molecule with the temperature and pressure dependent chemical potential, given as

where μ_O2(T, p) is the chemical potential of molecular oxygen at temperatures T and pressure p and μ_O2(T₀, p₀) is the total energy obtained from DFT calculations (at T = 0 K and p = 0 bar). At a given pressure, p₁, the temperature dependent contribution to the chemical potential is given as δ_O2(T, p₁) and is obtained from the JANAF thermochemistry tables,³⁸ and the last term is the pressure dependent contribution under the ideal gas assumption. Vibrational contributions to the entropy are not considered. The candidates for oxygen vacancies are calculated using the same method as for the desodiated phases, and only the atomic positions are allowed to relax in order to simulate vacancies in the dilute limit. To avoid self-interaction errors, the supercells are adjusted as needed to keep each lattice parameter above 10 Å.

The Atomic Simulation Environment (ASE) has been used for automisation of calculations and for data analysis.³⁹ All crystal structures are rendered with VESTA.⁴⁰

2 Structural stability

We performed a search for experimentally known and theoretically predicted compounds with the stoichiometry A₅MO₄ (A = alkali metals, M = p-, d- and f-block elements) in the ICSD and Materials Project databases to get a set of starting geometries for our investigation. The list of all starting geometries used are listed in Table S1 (ESI†).

For all starting geometries, the reference structures were subject to a simple ion substitution and volume change to obtain input structures with stoichiometry Na₅FeO₄ and the starting volume 147 ų f.u.⁻¹. The volumes were changed preserving the ratio between the lattice axes. Energy-volume curves were computed for each of the starting geometries after being pre-relaxed as described in Section 1. The equilibrium energies, volumes, bulk moduli and its pressure derivative for all structures are summarised in Table S2.† The relative equilibrium energies along with the energy–volume curves of the ten structures with the lowest energies are shown in Fig. 1.

Fig. 1 (a) Equilibrium energies of the ten lowest energy structures, with the three modifications chosen for further investigation highlighted in a darker colour. (b) Energy-volume curves fitted to the Birch–Murnagahn equation of state for the ten lowest energy trial structures.

The most stable phase was found to have Pmmn space group symmetry (s.g. 59), with two distinct structures with Pbca space group symmetry (s.g. 61) lying within 100 meV f.u.⁻¹ of the most stable phase. This energy difference was chosen as a cut-off for further investigations. In order of increasing energies, they are hereby designated as β-, α- and γ-Na₅FeO₄, respectively, the two former following the naming convention by Hoppe and co-workers^41,42 and the latter given by the current authors. γ-Na₅FeO₄, which is predicted as the least stable of the three polymorphs at 0 K, corresponds to the experimentally reported structure of Na₅FeO₄ reported by Brachtel and Hoppe in 1978 (ref. 13) and by Fruchart et al. in 1985,⁴³ and is shown later in this work to stabilise at higher temperatures when also including the entropic contributions. The α-phase corresponds to the experimentally reported structure of Li₅FeO₄ by Luge and Hoppe in 1984,⁴⁴ and the β-phase corresponds to the reported structures of Na₅InO₄ and Li₅AlO₄, also reported by Hoppe and co-workers.^41,45

The structural data for these three polymorphs are summarised in Table S4,† and the primitive unit cells are shown in Fig. S1.† All of these three phases are related to the antifluorite structure, and can be described as defect and distorted variants of the ideal structure, as shown in Fig. 2. The primitive unit cell of the α-phase corresponds to a 2 × 2 × 2 supercell of the antifluorite structure, while the primitive unit cell of β-Na₅FeO₄ can be transformed to arrive at the same antifluorite-type supercell, by employing a new set of lattice vectors: , and . The γ-phase can be thought of as a distorted 2 × 1 × 4 supercell of the ideal antifluorite structure, where only the top half of the supercell is shown in Fig. 2d. The antifluorite structure, exemplified by Na₂O (Fm3̄m, s.g. 225), has a cubic closest packing of O²⁻, with Na⁺ in all tetrahedral holes, as seen in Fig. 2a. In Na₅FeO₄, Fe³⁺ replaces one Na⁺ in each of the octants and is charge compensated by the creation of two Na⁺-vacancies (v′_Na). The three structures differ from each other in the ordering of Fe³⁺ and v′_Na relative to the other octants, and the degree of deviation of the oxygen sublattice from the ideal cubic closest packing. The β-phase deviates the least from this ideal packing and retains a readily identifiable packing sequence. So does the α-phase, but here the oxygen sublattice is further distorted. The γ-phase deviates the strongest from the ideal antifluorite structure with alternating depressions in the cubic octants, highlighted in Fig. S2,† making it hard to identify the relation to the antifluorite structure at first glance.

Fig. 2 (a) A 2 × 2 × 2 supercell of Na₂O (Fm3̄m). (b) Supercell of the β-polymorph (Pmmn) that has undergone a unit cell transformation with new lattice vectors , and . (c) Unit cell of the α-polymorph (Pbca). (d) Partial view of a supercell of the γ-polymorph (Pbca), with and cut in half along the [001]-direction to include the eight face-centered cubes of the distorted O²⁻-sublattice.

The ordering of Fe³⁺ in both the β- and γ-polymorphs is such that entire sheets of only Na⁺ can be identified, where Na⁺ can be thought to diffuse more freely, whereas in the α-phase, no such sheets can be isolated. In the case of the β-phase, these sheets are found parallel to the (100)- and (010)-planes, and in the γ-phase such sheets exist parellel to the (100)-planes.

In the ideal cubic antifluorite structure, all O–O-distances are the same, but as Fe³⁺-cations and Na⁺-vacancies are introduced and the oxygen sublattice distorted, these distances will start to differ depending on the relative proximity of each O²⁻ to Fe³⁺, Na⁺ or v′_Na, respectively. In the Li-analogue, Li₅FeO₄, all FeO₄-tetrahedra are still connected to neighbouring LiO₄-tetrahedra or v′_LiO4-tetrahedra, meaning that there are shorter O–O-distances connecting these tetrahedra than some of the internal O–O-distances within the FeO₄-tetrahedron itself. In Na₅FeO₄, however, the FeO₄-tetrahedra become completely isolated. For this reason, this structure type is also referred to as an island structure by Hoppe¹³ and denoted Na₅[FeO₄]. This is the case for all polymorphs considered here.

We further investigate the mechanical and dynamical stability of each of the polymorphs. First, in order to assess the mechanical stability, we computed the single crystal elastic constants. The resulting matrices were diagonalised and resulted in only positive valued eigenvalues, indicating that all considered polymorphs are mechanically stable. The eigenvalues are given in Table 1 while the full symmetrised elastic moduli are given in Table S3.† Phonon calculations were then carried out in order to assess the dynamical stability of each phase, and to investigate any changes in relative stability with temperature. The atom-projected phonon densities of states (PDOS) are given in Fig. 3, while the full dispersion relations are given in Fig. S3–S5.† As can be seen, there is an absence of any imaginary phonon modes in all polymorphs, indicating that they are all predicted to be dynamically stable. As the crystal structure of β-Na₅FeO₄ is less distorted than for the other phases, this phase is also expected to have more degenerate vibrational modes and thus more defined peaks, which is also clearly visible. In all cases, Fe and Na dominate at lower frequencies, but while contributions from Fe diminishes from 5–6 THz and upwards, the Na-contributions remain considerable for several THz above this. In this region in particular, the independence of Na from the other atoms is visible, indicative of decent diffusive properties, which are investigated further with NEB-calculations in Section 3. Above the phononic band gap there are many degenerate states composed of Fe and O, indicating stiff bonds of more covalent nature.

Table 1 Eigenvalues (GPa) resulting from the diagonalisation of the elastic moduli, showing only positive values indicating mechanical stability of all three phases

Phase	1	2	3	4	5	6
β-NFO	161.9	52.1	67.1	19.3	33.9	14.9
α-NFO	158.1	39.3	50.3	30.8	31.6	27.9
γ-NFO	144.6	44.7	59.4	27.6	28.7	26.1

---

Fig. 3 The site-projected phonon density of states (PDOS) normalised per atom for (a) the β-polymorph, (b) the α-polymorph and (c) the γ-polymorph.

In addition, the Helmholtz free energies were calculated from the phonon calculations as a function of temperature, given by F(T) = E_DFT + E_vib − TS, where E_DFT is the ground-state energy obtained through conventional energy minimisation and E_vib − TS are the vibrational contributions to the Helmholtz free energy. The results are shown in Fig. 4. We see that the experimentally reported γ-phase, although disfavoured at 0 K, has a larger energy gain with temperature and becomes stabilised with respect to the other polymorphs. Our calculations predict that γ-Na₅FeO₄ becomes the most stable polymorph at around 430 K, as can be seen in the inset of Fig. 4a and in Fig. 4b. From Fig. 1b, it can be seen that also at lower volumes (i.e. higher pressures), there is a shift in relative stabilities, where both the β- and α-phases are destabilised with respect to the γ-phase.

Fig. 4 (a) The temperature adjusted relative energies of the α-, β- and γ-polymorphs. Inset shows the crossover point where γ-Na₅FeO₄ becomes the most stable polymorph. (b) The difference in total energies as a function of temperature with the total energy of the β-phase set as the reference.

This compound has previously been synthesised via a high-temperature solid state reaction route,^13,43 but since these structures are so close in total energy, it could mean that the other polymorphs would be available through careful tuning of the synthesis parameters, as a function of particle size or through slight chemical modifications.

3 Ionic transport

The diffusivity of Na⁺-ions in the three most stable phases through vacancy diffusion was further calculated using the NEB-method as detailed in Section 1. Jumps to interstitial sites are for practical reasons not considered here. All atomic labels herein refer to the Wyckoff positions listed in Table S4,† and the directions refer to the antifluorite-type unit cells described in Section 2. A summary of the diffusion paths identified, as well as a summary of each individual jump, can be found in Tables S5 and S6.†

For the β-phase, four unique Na–Na-jumps were identified giving rise to five unique diffusion pathways, two having a net movement in the [110]-direction, one along [001] and two with a net movement along [11̄0]. The most facile path has an overall barrier of 0.207 eV with a net movement in the [110]-direction involving two distinct jumps between Na1. Based on the limiting energy barrier along this path (0.207 eV), and the Na–Na-distance associated with this jump (2.79 Å), the diffusion coefficient along this direction is calculated according to eqn (4) to be 2.66 × 10⁻⁶ cm² s⁻¹ at 300 K, being by far the highest diffusion coefficient of all the polymorphs. The diffusion pathway and associated energy barriers are shown in Fig. 5a–c. The movement is contained within a sheet free of any Fe³⁺-ions that exist in this phase due to the particular ordering of the transition metal cations. The remaining paths have barriers that range from 0.301–0.406 eV. Of these, two are also purely contained within sheets of Na, while the two others have a more tortuous path that involves movement between planes.

Fig. 5 (a) The most facile diffusion path through the antifluorite-type unit cell, (b) a cutout of the diffusion path and (c) the associated migration barriers calculated from NEB-calculations, respectively, of the β-polymorph. (d–f) shows the same for the α-polymorph. All labels refer to the Wyckoff site occupation tabulated in Table S4.†

For the α-polymorph, ten unique jumps were identified which give rise to four percolating pathways, with two having a net movement in the [001]-direction, one along [010] and the last along [100]. Here, the lowest overall barrier was calculated to be 0.434 eV along [001], following the pathway Na2–Na4–Na3, shown in Fig. 5d–f. The diffusion coefficient is calculated to be 1.4 × 10⁻⁹ cm² s⁻¹ based on the limiting jump along this pathway, between Na2 and Na4, with a jump distance of 2.88 Å and an energy barrier of 0.425 eV. This is, notably, three orders of magnitude lower than what was found in the β-phase. The diffusion path along the [010]-direction, and the other path along [001], are limited by the same jump, and have calculated diffusion coefficients of 1.08 × 10⁻¹⁰ cm² s⁻¹. The last diffusion path along [100] involved a jump that is not well-described as a jump between two Na-sites, and as diffusion via interstitials and/or interstitialcies is not treated here, the reported values are based on the jump with the highest barrier under the assumption that this truly is the limiting jump along the diffusion path.

For the γ-phase, four unique jumps were identified giving rise to two unique pathways along the [100]- and [001]-directions, respectively. However, none of these pathways consisted purely of jumps from one Na-site to another, which indicates that diffusion through interstitials and/or interstitialcies must be considered for an accurate description of diffusive properties of γ-Na₅FeO₄. However, along both pathways, there is a relatively large barrier (0.364 eV) that is assumed to be the limiting jump in order to give an estimation of the diffusive properties in this polymorph.

In Fig. S6† the diffusion coefficients for the most facile diffusion paths for each polymorph are shown as a function of temperature as calculated from eqn (4). The values of the diffusion coefficients are quite high compared to experimental values reported for other known SIB cathode materials.^46–49 If the values predicted in this work could be achieved experimentally, especially if the β-phase could be realised, Na₅FeO₄ could show excellent rate capabilities suitable for high power applications.

4 Structural stability upon desodiation

In order to assess a likely path for stable electrochemical cycling, we compared desodiated structures based on the three pristine polymorphs with other structure types obtained through a similar search as in Section 2 for Na_5−xFeO₄ (x = 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0). As the number of desodiated configurations scales with the binomial distribution as x increases, with a peak at x = 2.5, only a subset of all possible configurations was considered as outlined in Section 1. Since O₂-outgassing might be expected in Na-rich cathode materials during charging, we also considered oxygen deficient structures given by Na_5−xFeO_4−y (y = 0.5, 1.0, 1.5, 2.0). These structure types were also obtained through a similar structure search as above. We refer to desodiated variants derived from the pristine phases by their pristine label, while all other structure types are grouped together based on oxygen stoichiometry.

For each level of desodiation, x, we compare the total energies between these groups to evaluate the reversibility of desodiation, where the total energies of the oxygen deficient modifications are calculated from eqn (2). We considered a desodiation step irreversible once oxygen deficient structures are preferred. The modifications with the lowest total energies from each group are plotted in Fig. 6, and it can be seen that up to 1.5 Na⁺ can be extracted before phases with lower oxygen content are favoured. Based on this, we found that there are two intervals of stable cycling without associated oxygen loss: (1) Na₅FeO₄ ⇌ Na_3.5FeO₄ + 1.5Na⁺ and (2) Na_3.5FeO₃ ⇌ Na₂FeO₃ + 1.5Na⁺. To access the second interval, there must necessarily be a formation step with initial oxygen loss to arrive at Na_3.5FeO₃ before the subsequent cycling can occur without further loss of oxygen. Such a cycling window could be of interest, where the initial cycling window could serve as a presodiation step, although the associated O₂-evolution might pose safety hazards.

Fig. 6 Relative energies of Na_5−xFeO_4−y (x = 0.5–3.0, y = 0.5–2.0). Each panel corresponds to a given level of desodiation, and each bar corresponds to the lowest energy modification within its respective group of trial structures. Three of the groups are desodiated variants of the β-, α- and γ-polymorphs, respectively, while the rest are made up of structures from a database search for the given stoichiometry. Where there are no bars, no structure matching the stoichiometry showed up in the search. For each level of desodiation, the energies are given relative to the lowest energy modification.

Convex hulls were constructed for these two intervals according to eqn (1), and are shown in Fig. 7. For the first interval we identified two plateaus corresponding to Na₅FeO₄ → Na₄FeO₄ + 1Na⁺ and Na₄FeO₄ → Na_3.5FeO₄ + 0.5Na⁺. The associated structures of Na₄FeO₄ and Na_3.5FeO₄ have triclinic (P1̄, s.g. 2) and monoclinic (P2₁, s.g. 4) space group symmetries, respectively, and structural data for these structures is tabulated in Table S7.† Both phases are still related to the pristine phase, preservering the isolated FeO₄-tetrahedra, and retaining the general packing of the oxygen sublattice, although with further distortions. This indicates that this cycling interval should not be prohibited by excessive structural rearrangements. The associated volume changes along this desodiation path is shown in Fig. 8a, showing a volume contraction of just under 18%. This contraction is on the same order of magnitude as the volume expansion in the common SIB anode material hard carbon,⁵⁰ providing a good match to reduce the overall stress associated with these volume changes in the cell during cycling. The predicted equilibrium voltages for these plateus were calculated using the phases along the convex hull according to eqn (3), and found to be 1.78 V for the first plateu and 2.76 V for the second, corresponding to an average voltage of 2.11 V for the whole interval, shown in Fig. 9a. With the exchange of 1.5 e⁻, this gives a specific capacity of 171.23 mA h g⁻¹ and a gravimetric energy density of 361.30 W h kg⁻¹.

Fig. 7 Convex hull constructred for the reactions (a) Na₅FeO₄ ⇌ Na_3.5FeO₄ + 1.5Na⁺ and (b) Na_3.5FeO₃ ⇌ Na₂FeO₃ + 1.5Na⁺.

Fig. 8 Relative volume changes corresponding to the reactions (a) Na₅FeO₄ ⇌ Na_3.5FeO₄ + 1.5Na⁺ and (b) Na_3.5FeO₃ ⇌ Na₂FeO₃ + 1.5Na⁺. Volume changes of desodiated variants of the β-, α- and γ-polymorphs are plotted for comparison in faded colours. The volume changes are given with respect to the initial volume for the given series.

Fig. 9 Predicted potential profiles vs. Na⁺/Na for the reactions (a) Na₅FeO₄ ⇌ Na_3.5FeO₄ + 1.5Na⁺ and (b) Na_3.5FeO₃ ⇌ Na₂FeO₃ + 1.5Na⁺.

For the second interval, a convex hull was similarily constructed (Fig. 7b) including only oxygen deficient phases, to reveal another two plateaus: Na_3.5FeO₃ → Na₃FeO₃ + 0.5Na⁺ and Na₃FeO₃ → Na₂FeO³ + 1Na⁺. The corresponding potential profile (Fig. 9b) exhibits a very low voltage of 0.15 V for the first plateu and 2.07 V for the second. The very low voltage of the first step limits the practical use to the interval Na₃FeO₃ ⇌ Na₂FeO₃ + 1Na⁺. This corresponds to a specific capacity of 114.15 mA h g⁻¹ and an energy density of 236.29 W h kg⁻¹ for the limited cycling window. In this case, the initial phase (Na_3.5FeO₃) has triclinic (P1, s.g. 1) space group symmetry, while the desodiated phases (Na₃FeO₃ and Na₂FeO₃) both have monoclinic (P2₁/c, s.g. 14 and C2/m, s.g. 12) space group symmetries, for which structural data is tabulated in Table S8.† These structures differ considerably from each other, and the Fe-coordination changes from edge-sharing tetrahedra, through corner-sharing tetrahedra, to edge-sharing octahedra in the last stage. This would require quite large reconstructions between each de-/sodiation step and is likely kinetically prohibitive. In addition, the volume changes associated with this cycling interval (Fig. 8b) are much larger than for the first interval. For these reasons, this cycling interval is not further investigated.

To look further into the stability of desodiation with respect to oxygen loss, we calculated the oxygen vacancy formation enthalpies for different levels of desodiation along the convex hull and for each of the three phases, according to the method outlined in Section 1. The most stable configurations of the desodiated variants previously obtained were used as the basis for these calculations. The vacancy concentrations for each considered modification, along with the vacancy formation enthalpies at 0 K, are summarised in Table S9.†

In Fig. 10 the evolution of the vacancy formation energies as a function of desodiation level x for the phases found along the convex hull is given for p_O2 = 0 bar at T = 0 K (Fig. 10a) and adjusted for ambient conditions, p_O2 = 0.2 bar and T = 300 K (Fig. 10b). The evolution of the vacancy formation energy with temperature is given for each individual composition in Fig. S7.† Oxygen vacancies are heavily unfavoured in the fully discharged state. However, as the structures are desodiated, the vacancies are quickly stabilised. Along the convex hull, the formation of an oxygen vacancy is barely unfavoured at 0 K at x = 1.5 (0.05 eV), however adjusting for temperature and pressure, the vacancy becomes slightly favoured at this level of desodiation (−0.28 eV). This is consistent with what is shown in Fig. 6, where we saw that oxygen deficient phases become stable after this level, and these results show that we can begin to expect a small amount of irreversible oxygen loss already at x = 1.5. As the vacancy concentrations are quite low, this should not be prohibitive, as any further oxygen loss at this point is found to be unfavourable. Furthermore, this decrease in vacancy formation energies with desodiation is also an indication that a certain amount of charge compensation can be attributed to oxygen, and that there is reversible contribution of anionic redox up until a point where oxide ions are fully oxidised and irreversibly disassociates from the lattice. The possible charge compensation mechanism will be expanded upon in the next section.

Fig. 10 The evolution of vacancy formation energies as a function of x in Na_5−xFeO₄ at (a) p_O2 = 0 bar and T = 0 K, and (b) p_O2 = 0.2 bar and T = 300 K.

As expected, for the second cycling plateau, for which the compound has already suffered considerable oxygen loss, further vacancy formation is again heavily unfavoured at the start of the cycling interval. It also trends towards stabilisation upon further desodiation, however, it is never predicted to become stable at ambient conditions within the cycling interval.

5 Electronic structure and charge compensation

We herein report on the electronic structure of the pristine compound and the phases along the convex hull in the cycling interval Na₅FeO₄ → Na₄FeO₄ + 1Na⁺, and further discuss the charge compensation mechanism upon desodiation in this interval. In Fig. 11, the site and orbital projected electronic density of states (OPDOS) is given for each species for the pristine compound Na₅FeO₄, along with the crystal orbital overlap population (COOP) of the Fe–O and O–O interactions within the FeO₄-tetrahedra. By convention, the Fermi-level (E_F) is set to 0 eV. The expected features from a molecular orbital approach of a tetrahedrally coordinated Fe³⁺ compound with a high on-site repulsion of the d⁵-configuration are reproduced. Below E_F, the t₂-and t₂*-states show a mixed Fe/O-character, although the O-states dominate, making this a charge transfer band gap. The band gap is computed to be just below 2 eV, however the number of available states at the bottom of the conduction band is quite small, and the majority of available states are found in the empty Fe 3d-states at about 3.5 eV above E_F. This can be more clearly seen in the full bandstructure, given in Fig. S8.† This should be non-prohibitive from an electrochemical perspective, as shown by the commercial success of LiFePO₄ that has a band gap exceeding this,⁵¹ where the combination of conductive additives and particle size reductions are used as means to facilitate electron transport to the redox centres,⁵² and a similar strategy could be pursued for Na₅FeO₄.

Fig. 11 The site and orbital projected electronic density of states (OPDOS) normalised per atom for Na, Fe and O, respectively, along with the crystal orbital overlap population (COOP) computed for Fe–O and O–O interactions within the FeO₄-tetrahedra in β-Na₅FeO₄. The spin up channel is plotted with positive values and the spin down channel is plotted with negative values in the OPDOS-plots. Bonding interactions are plotted with positive values and antibonding interactions are plotted with negative values in the COOP-plots. The Fermi-level (E_F) is by convention set to 0 and is indicated by a vertical, dashed black line.

The low contribution of Na states around E_F confirms the ionic nature of the Na–O interactions, while the high degree of degeneracy of Fe 3d- and O 2p-states indicate a rather large degree of covalency in the Fe–O interactions, as previously seen in the phononic density of states. As is the case in other Li/Na-rich transition metal oxides, the environments surrounding each oxygen lead to localised O 2p states, often referred to as O 2p “lone pairs”. The oxide ions in Na₅FeO₄ are not cubically coordinated as in the ideal antifluorite structure, but instead octahedrally coordinated, with the octahedra being heavily distorted. In these octahedra, due to the presence of Na–O–Na bridges, O 2p-orbitals pointing along these directions are left essentially unhybridised as lone pairs on the oxide ions. The presence of these lone pairs can be identified through the computation of the electron localisation function (ELF), a topological tool of the electron density that identifies regions in space where paired electrons are localised.⁵³ Indeed, electron pairs can be shown to exist along these directions as shown in Fig. S13a.† These localised oxygen states are also evident in the OPDOS, where one pair can be found as sharp peaks in both spin channels at around 1 eV below E_F and the other at around 2 eV below, although here not as readily identifiable, as these are not quite as localised and overlaps the t₂-states. Due to the rather large distortion of the structure, there is bound to be a certain degree of overlap that would otherwise be symmetry disallowed, and this could lead to a lower degree of localisation for the lone pair states at lower energies. Interestingly, these lone pairs that are otherwise widely accepted to be the origin of anionic redox in other Li/Na-rich compounds,⁸ seems to be inactive upon desodiation, as seen in Fig. S13b and c.† The localised electron pairs are still present even at x = 1.5, with no apparent change at the same isosurface values.

This perserverance of the highly localised O 2p-states are also evident in the OPDOS of the desodiated phases, as seen in Fig. 12, however, as the structure becomes more distorted, the peaks become less sharply defined. Instead, “hole” states composed of a combination of empty Fe 3d- and O 2p-states become evident above E_F, first one at x = 1.0 and then a second one at x = 1.5. The Fe 3d-states make up about a quarter of the total number of states in each of these holes. Upon integration of the hole states at x = 1.0, we find that the hole corresponds to a single electron per formula unit, as expected. For x = 1.5, the hole closest to E_F integrates to one electron per formula unit and the hole at a slightly higher energy level sums up to half an electron per formula unit. In such a situation we would expect a charge disproportionation where half the redox centres will have donated two electrons for charge compensation, while the other half will have donated only one. This is also what is observed from visualisation of the electronic states associated with each hole at different FeO₄-tetrahedra, as seen in the inset of Fig. 12. This is also evident from inspection of the structural changes, shown in Fig. S14:† for x = 0.0 and 1.0 there is only a single distinct Fe-site, while for x = 1.5 there are four. Among those four, two and two and pairwise nearly identical and are subsequently labelled Fe1/Fe2 and Fe3/Fe4. One of these pairs is very similar to the FeO₄-tetrahedron found in x = 1.0, while the other is further changed. The Fe–O distances are lowered for each desodiation step, from 1.93–1.94 Å at x = 0.0, to 1.81–1.83 Å at x = 1.0, to 1.73–1.77 Å for the most oxidised tetrahedra at x = 1.5. These values match well with what is expected from using ionic radii of Shannon⁵⁴ of Fe with oxidation states +3, +4 and +5 (obtained by interpolating values for tetrahedrally coordinated Fe between Fe³⁺ and Fe⁶⁺): 1.89 Å, 1.81 Å and 1.73 Å. The distances between oxygen ions within each tetrahedron also decreases with desodiation, and this is further discussed below.

Fig. 12 The site and orbital projected electronic density of states (OPDOS) for O 2p and Fe 3d for Na₄FeO₄ (x = 1.0) and Na_3.5FeO₄ (x = 1.5). Insets show the partial charge densities of the symmetry equivalent FeO₄-tetrahedra associated with the “hole” states that appear upon desodiation. The labels of Fe corresponds to Wyckoff positions listed in Table S4.† For Na_3.5FeO₄, Fe1/Fe2 and Fe3/Fe4 are pairwise almost identical, and treated as two distinct Fe-sites.

To further investigate the nature of the redox in Na₅FeO₄, we computed the Fukui function at all desodiation levels. To compute these, we keep the structures frozen, and recompute the charge densities with either one electron added (Fukui+) or subtracted (Fukui−) from the system, revealing the spatial localisation of the lowest unoccupied and the highest occupied electronic levels respectively. A visual representation of the results from these calculations is shown in Fig. 13. This reveals that along the whole cycling interval, the electronic levels associated with the charge compensation are hybridised Fe/O-states, indicating that all redox activity is contained to these levels, and again shows that the O 2p lone pairs are left intact.

Fig. 13 The Fukui functions computed for the phases along the convex hull. This shows the electronic reorganisation upon addition (Fukui+) and removal (Fukui−) of one electron to the system, being visible as a decrease (blue) or an increase (green) of electronic density. The isosurface value is 0.002 in all cases. Atomic labels refer to Wyckoff positions given in Table S4.† Fe1/Fe2 and Fe3/Fe4 are pairwise almost identical to each other, and are treated as two distinct Fe-sites.

The emptying of electron states from the covalent Fe–O bonds is also predicted to occur in the Li-rich cathode material Li_1.2Nb_0.3Fe_0.4O₂ with the simultaneous formation of a superoxide species (O⁻₂) between two oxide ions of neighbouring octahedra, with a bond length of 1.36 Å.⁵⁵ This is not a stable situation, and leads to a reductive coupling with the oxidised Fe-cation, and leads to irreversible oxygen loss. While shortening of O–O distances is also predicted upon desodiation of Na₅FeO₄, these distances are still far above what could be considered superoxide or even peroxide species (O²⁻₂), with the shortest O–O separation at x = 1.5 being 2.67 Å, as seen in Fig. S14.† This is more in line with what is observed in another tridimensional model system, β-Li₂IrO₃ (2.53 Å), that exhibits reversible electrochemical cycling without any oxygen loss.^56,57

6 Conclusion

The sodium-rich oxoferrate Na₅FeO₄ was investigated as a potential cathode for SIBs from first principles using density functional theory. Out of multiple starting geometries, we find three low-energy orthorhombic polymorphs related to the antifluorite structure lying within 100 meV f.u.⁻¹ of each other, including the experimentally reported structure that is found to be favoured at temperatures below typical synthesis temperatures upon inclusion of entropic contributions. All polymorphs are found to be both mechanically and dynamically stable. NEB calculations point toward excellent diffusive properties, although further work is needed to investigate the role of interstitial and/or interstitialcy diffusion in these materials. Investigations are made into relative phase stability as sodium is removed from the structure. We have looked into direct desodiation of the pristine structures, phase transformations and possible oxygen loss, revealing the existence of two possible reversible cycling intervals: Na₅FeO₄ ⇌ Na_3.5FeO₄ + 1.5Na⁺ and Na_3.5FeO₃ ⇌ Na₂FeO₃ + 1.5Na⁺. While the first interval shows promise with a specific energy density of more than 360 W h kg⁻¹, the second interval is limited by too severe reconstructions to be a likely candidate as a cathode material. Analysis of the electronic structure, including the computation of the electron localisation function and the Fukui functions, point towards a redox mechanism contained to hybridised Fe 3d-/O 2p-levels where FeO₄-tetrahedra serve as a joint redox centre during de-/sodiation, while unhybridised O 2p “lone pairs”, usually reported to be involved in the anionic redox process, are left untouched.

If realised experimentally, Na₅FeO₄ has the potential of becoming a very inexpensive cathode material for use in stationary storage applications, consisting of only cheap and geographically abundant raw materials. The results presented show that it warrants further experimental investigations, and such investigations are currently underway.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors gratefully acknowledge the Research Council of Norway for providing the computer time (under the project numbers NN2875k and NS2875k) at the Norwegian supercomputer facility.

References

1. C. Delmas, Adv. Energy Mater., 2018, 8, 1–9.
2. European Commission, Critical Raw Materials Resilience: Charting a Path towards greater Security and Sustainability, European commission technical report, 2020.
3. C. Zu, Y. Ren, F. Guo, H. Yu and H. Li, Adv. Energy Sustain. Res., 2021, 2, 2100062.
4. B. Swain, Sep. Purif. Technol., 2017, 172, 388–403.
5. A. Alessia, B. Alessandro, V. G. Maria, V. A. Carlos and B. Francesca, J. Cleaner Prod., 2021, 300, 126954.
6. P. Xu, D. H. Tan and Z. Chen, Trends Chem., 2021, 3, 620–630.
7. G. Assat and J. M. Tarascon, Nat. Energy, 2018, 3, 373–386.
8. M. Ben Yahia, J. Vergnet, M. Saubanère and M. L. Doublet, Nat. Mater., 2019, 18, 496–502.
9. M. Sathiya, G. Rousse, K. Ramesha, C. P. Laisa, H. Vezin, M. T. Sougrati, M. L. Doublet, D. Foix, D. Gonbeau, W. Walker, A. S. Prakash, M. Ben Hassine, L. Dupont and J. M. Tarascon, Nat. Mater., 2013, 12, 827–835.
10. S. Narukawa, Y. Takeda, M. Nishijima, N. Imanishi, O. Yamamoto and M. Tabuchi, Solid State Ionics, 1999, 122, 59–64.
11. A. Hirano, T. Matsumura, M. Ueda, N. Imanishi, Y. Takeda and M. Tabuchi, Solid State Ionics, 2005, 176, 2777–2782.
12. C. Zhan, Z. Yao, J. Lu, L. Ma, V. A. Maroni, L. Li, E. Lee, E. E. Alp, T. Wu, J. Wen, Y. Ren, C. Johnson, M. M. Thackeray, M. K. Chan, C. Wolverton and K. Amine, Nat. Energy, 2017, 2, 963–971.
13. G. Brachtel and R. Hoppe, Z. Anorg. Allg. Chem., 1978, 446, 77–86.
14. W. Burow and R. Hoppe, Naturwissenschaften, 1980, 67, 192.
15. H. Zentgraf and R. Hoppe, Z. Anorg. Allg. Chem., 1980, 462, 61–70.
16. G. Brachtel, N. Bukovec and R. Hoppe, Z. Anorg. Allg. Chem., 1984, 515, 101–113.
17. A. Belkly, M. Helderman, V. L. Karen and P. Ulkch, Acta Crystallogr., Sect. B: Struct. Sci., 2002, 58, 364–369.
18. A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder and K. A. Persson, APL Mater., 2013, 1, 011002.
19. P. E. Blöchl, Phys. Rev. B, 1994, 50, 17953–17979.
20. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
21. G. Kresse and J. Hafner, Phys. Rev. B, 1993, 47, 558–561.
22. G. Kresse and J. Hafner, Phys. Rev. B, 1994, 49, 14251–14269.
23. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
24. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
25. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
26. S. Dudarev and G. Botton, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 57, 1505–1509.
27. P. E. Blöchl, O. Jepsen and O. K. Andersen, Phys. Rev. B, 1994, 49, 16223–16233.
28. M. Hebbache and M. Zemzemi, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 70, 5–10.
29. A. Togo and I. Tanaka, Scr. Mater., 2015, 108, 1–5.
30. C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke and T. E. Oliphant, Nature, 2020, 585, 357–362.
31. A. Urban, D. H. Seo and G. Ceder, npj Comput. Mater., 2016, 2, 16002.
32. S. P. Ong, W. D. Richards, A. Jain, G. Hautier, M. Kocher, S. Cholia, D. Gunter, V. L. Chevrier, K. A. Persson and G. Ceder, Comput. Mater. Sci., 2013, 68, 314–319.
33. L. Ward, A. Dunn, A. Faghaninia, N. E. Zimmermann, S. Bajaj, Q. Wang, J. Montoya, J. Chen, K. Bystrom, M. Dylla, K. Chard, M. Asta, K. A. Persson, G. J. Snyder, I. Foster and A. Jain, Comput. Mater. Sci., 2018, 152, 60–69.
34. G. Henkelman and H. Jónsson, J. Chem. Phys., 2000, 113, 9978–9985.
35. G. Henkelman, B. P. Uberuaga and H. Jónsson, J. Chem. Phys., 2000, 113, 9901–9904.
36. E. Bitzek, P. Koskinen, F. Gähler, M. Moseler and P. Gumbsch, Phys. Rev. Lett., 2006, 97, 1–4.
37. K. Kang, D. Morgan and G. Ceder, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 1–4.
38. M. Chase, JANAF Thermochemical Tables, American Chemical Society, American Institute of Physics for the National Bureau of Standards, New York, NY, 3rd edn, 1986.
39. A. Hjorth Larsen, J. Jørgen Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. Bjerre Jensen, J. Kermode, J. R. Kitchin, E. Leonhard Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. Bergmann Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiøtz, O. Schütt, M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. Zeng and K. W. Jacobsen, J. Phys.: Condens. Matter, 2017, 29, 273002.
40. K. Momma and F. Izumi, J. Appl. Crystallogr., 2011, 44, 1272–1276.
41. H. Glaum and R. Hoppe, Z. Anorg. Allg. Chem., 1990, 583, 24–30.
42. F. Stewner and R. Hoppe, Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem., 1971, 27, 616–621.
43. D. Fruchart, J. L. Soubeyroux, H. Kessler and J. M. Lassalle, J. Solid State Chem., 1985, 57, 191–196.
44. R. Luge and R. Hoppe, Z. Anorg. Allg. Chem., 1984, 513, 141–150.
45. R. Hoppe and H. König, Z. Anorg. Allg. Chem., 1977, 430, 211–217.
46. Y. Qi, F. Xiao, H. Wang, Y. Yan, S. Bao and M. Xu, Chem. Eng. J., 2021, 405, 127035.
47. M. H. Han, E. Gonzalo, G. Singh and T. Rojo, Energy Environ. Sci., 2015, 8, 81–102.
48. D. Tie, G. Gao, F. Xia, R. Yue, Q. Wang, R. Qi, B. Wang and Y. Zhao, ACS Appl. Mater. Interfaces, 2019, 11, 6978–6985.
49. M. Chen, Y. Zhang, G. Xing and Y. Tang, Front. Chem., 2020, 8, 1–7.
50. A. Agrawal, S. Janakiraman, K. Biswas, A. Venimadhav, S. K. Srivastava and S. Ghosh, Electrochim. Acta, 2019, 317, 164–172.
51. F. Zhou, K. Kang, T. Maxisch, G. Ceder and D. Morgan, Solid State Commun., 2004, 132, 181–186.
52. S. Ferrari, R. L. Lavall, D. Capsoni, E. Quartarone, A. Magistris, P. Mustarelli and P. Canton, J. Phys. Chem. C, 2010, 114, 12598–12603.
53. A. D. Becke and K. E. Edgecombe, J. Chem. Phys., 1990, 92, 5397–5403.
54. R. D. Shannon and C. T. Prewitt, Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem., 1970, 26, 1046–1048.
55. N. Yabuuchi, M. Nakayama, M. Takeuchi, S. Komaba, Y. Hashimoto, T. Mukai, H. Shiiba, K. Sato, Y. Kobayashi, A. Nakao, M. Yonemura, K. Yamanaka, K. Mitsuhara and T. Ohta, Nat. Commun., 2016, 7, 1–10.
56. E. McCalla, A. M. Abakumov, M. Saubanère, D. Foix, E. J. Berg, G. Rousse, M. L. Doublet, D. Gonbeau, P. Novák, G. Van Tendeloo, R. Dominko and J. M. Tarascon, Science, 2015, 350, 1516–1521.
57. P. E. Pearce, A. J. Perez, G. Rousse, M. Saubanère, D. Batuk, D. Foix, E. McCalla, A. M. Abakumov, G. Van Tendeloo, M. L. Doublet and J. M. Tarascon, Nat. Mater., 2017, 16, 580–586.

---

Footnotes

† Electronic supplementary information (ESI) available. See https://doi.org/10.1039/d2ra01834a

‡ Present address: Centre for Bioinformatics, University of Oslo, P.O. Box 1080 Blindern, N-0316 Oslo, Norway.

---