Ab initio random structure searching for battery cathode materials

Lithium-ion batteries (LIBs) are electrochemical energy storage devices characterized by high energy density, good rate capability, and a long shelf-life. They have dominated the market for portable electronics since their commercialization in the 1990s.^1,2 This technology has been challenged by the increasing demands for higher energy density, due to the growing electric vehicle market, and grid-scale energy storage. One of the most important performance parameters, the energy density of a battery, is determined by its average discharge voltage, lithium storage capacity, and weight (or volume). Among these factors, both the voltage and the capacity depend heavily on the physiochemical properties of cathode materials used. Therefore, the search for cathode materials with high discharge voltages and large capacities is one of the central aims of battery research.^3,4 Beyond that, due to the demand for high-power output by electric vehicles and grid-scale energy storage, a high-rate capability of a battery is also essential. In this context, cathodes that transport both lithium and electrons rapidly are critically needed.

Conventionally, new cathode materials are discovered through ad hoc design based on the understanding of known chemistry and crystallography. For example, the intercalation reaction between guest ions with solid hosts has been known for over a century.⁴ By studying the fundamental properties of the Li_xTiS₂ series, Whittingham used the layered TiS₂ as a cathode and demonstrated the first rechargeable lithium battery.⁵ However, such cathodes suffer from low redox potential vs Li/Li⁺, so the energy density must be improved. To achieve this, Goodenough and co-workers replaced S²⁻ with O²⁻ and enabled the use of the Co^3+/4+ redox couple based on the knowledge that the S²⁻ 3p band lies at a higher energy as compared with that of O²⁻ 2p.⁶ This led directly to the discovery of LiCoO₂, one of the most successful cathode materials for LIBs to date.⁷ Despite the commercial success, the search for cathode materials based on ad hoc design and experimental trial-and-error has proven difficult. In fact, only a very limited number of materials have been confirmed to be competitive, despite the tremendous efforts made in searching for them. A new paradigm for cathode material discovery is critically needed.

Computational materials science has become a powerful tool in materials discovery.^8–10 Thanks to developments in electronic structure methods, in particular, density functional theory (DFT) and related computational approaches, the first-principles computation of the properties of practical materials has become a routine in materials research.^11–14 At the same time, the dramatic increase in computational capacity of our research infrastructure has empowered researchers to carry out high-throughput screening of materials for those with desired properties.⁸ Databases have been assembled to archive properties such as formation energies and bandgaps.^15,16 In the context of cathode materials, phase stability, discharge voltage, and cation diffusion properties can now be calculated and analyzed based on these precomputed data entries.^17,18 Pioneering efforts have also been made to predict new cathode materials. For example, elemental substitution^19,20 and motif templating^21–24 have been proposed as a means to effectively generate candidate structures, followed by property screening. While such approaches led to the discovery of a number of cathodes,^19,20,25 the data records typically are limited in number. The compositional and structural space explored is severely constrained, and many possible combinations of chemical species have been overlooked. More importantly, even for a single given composition, crucial information is omitted that would aid in the exploration of new materials, namely, polymorphism, or different low energy structures. These metastable materials can sometimes be synthesized and display desirable properties. In fact, materials with polymorphism have been widely studied in a number of areas including high pressure physics, organic chemistry, and electrochemistry.^26–28 In the context of battery cathode materials, such polymorphism may provide exciting opportunities for materials with extraordinary electron and ion transport characteristics.

Here, we address these issues by showing how ab initio random structure searching (AIRSS) can be used as an efficient tool for the exploration of novel cathode materials for batteries.^29–31 The method is based on a random sampling of the first principles potential energy surfaces (PESs) of atomic arrangements and has been successfully applied in a number of areas including high pressure physics, superconductors, semiconductors, and energy storage.^30,32–35 Here, we demonstrate that by constraining the search space through the choice of parameters such as chemically aware minimum interatomic distances, cell volumes, and symmetries, the efficient identification of both stable and metastable cathode materials is possible. Specifically, we explore existing cathode compositions such as LiCoO₂, LiFePO₄, and Li_xCu_yF_z to showcase this method and show how to achieve an efficient sampling of the PES. Building on this, we apply the AIRSS method to lithium-stuffed transition-metal oxalates in the search for novel cathode materials. Screening the low energy outcomes of the search, we identify several oxalate polymorphs as good candidates with decent rate capability and energy density.

The general searching framework is illustrated in Fig. 1. It consists of the choice of composition, the crystal structure search using AIRSS, and the post-search screening. Composition-wise, we chose LiCoO₂, LiFePO₄, Li_xCu_yF_z, and Li₂TM(C₂O₄)₂ (TM = Fe, Co, Ni, V, and Mn) in this study. LiCoO₂ and LiFePO₄ were selected because they are the most representative materials of two cathode families: simple transition metal oxides and polyanionic compounds. It is worth noting that, to determine the stability of a phase, not only the composition of interest needs to be considered but also all the other possible ones within the given chemical space. While for most systems, the stable phases can be directly taken from databases, for systems that are less studied, an extensive search over the entire compositional space is necessary to construct the relevant phase diagram. In this case, we demonstrate the search of a phase diagram using Li_xCu_yF_z because Cu₂F is a high-capacity conversion-type cathode, which undergoes interesting conversion reactions during charge and discharge. Finally, Li₂TM(C₂O₄)₂ was chosen to highlight the capability of the current searching scheme because Li₂Fe(C₂O₄) has been reported to be a potential cathode with interesting anion redox properties while all the other members are less well studied.

After choosing the composition, crystal structure predictions were carried out using AIRSS. The details of this approach are provided in Subsection II B, and the impact of selecting the few search parameters on its efficiency is discussed in Sec. III A. After the AIRSS search, candidate structures were selected for property screening based on their phase stability. For a cathode material, the properties of interest are the voltage, the capacity, the cation transport barrier, and the electronic bandgaps. The first two determine the energy density of a cathode, while the latter two provide a proxy to predict the rate capability. These property calculations are standard in battery research and are briefly discussed in Sec. II C.

At the core of the searching scheme, AIRSS samples the local energy minima of a PES by first generating random “sensible” structures and then locally optimizing their geometries using first-principles calculations. The efficiency or even the success of the algorithm, as we will explain later, is highly dependent on choosing sensible parameters and constraints based on physical quantities, e.g., reasonable density, no close contacts, and proper neighbor relations. The general workflow of an AIRSS search is depicted in Fig. 1(a). Generating random structures with few constraints will provide the widest coverage of the search space and therefore reduce the possibility of completely missing out the global minimum. However, it will also result in low search efficiency, i.e., the rate of encountering lower energy structures is reduced. Considering that cathode materials operate at ambient pressures, it is natural to use chemical ideas to constrain the search space and steer the search toward chemically “sensible” regions. This is done by applying constraints during structure generation. To generate a structure that is chemically sensible, both the bond lengths and the neighboring relation between atoms need to be constrained. Two parameters are effective in this context. One is the minimum separation between atoms (which may be different for each pair of species and so chemically aware). To acquire structures that fulfill the distance constraints, they are generated randomly and are optimized using a hard sphere potential to reach the desired species-dependent minimum separations. It is important to note that such minimum separations are species-pairwise and an example is shown in Table I. Such distance constraints avoid close contacts between atoms, which may lead to numerical problems for the first principles code during the geometry relaxation. Furthermore, by fixing the cell volume (per formula unit/f.u.), such species-pairwise minimum separations can be used to avoid unwanted neighboring relations. For example, by setting the minimum separations between Li and Li and between Li and O to 2.8 and 2.0 Å, respectively, in a lithium containing oxide, OLi_x clusters will be favored instead of Li_x neighboring. Should one have knowledge on the chemical system, the values of the minimum separations and cell volumes can be measured from known structures. For the case where the chemistry is new, one can perform a preliminary AIRSS search using random minimum separations and then measure the resulting values from the structure with the lowest energy. In practice, to avoid overly constraining the system, the structural rejection criteria based on the minimum separation matrix are slightly loosened. Symmetry is also exploited during the AIRSS search as most known crystals are characterized by a degree of symmetry. Moreover, adopting symmetry accelerates the DFT calculations through reducing the number of k-points required to sample the Brillouin zone and the number of degrees of freedom in the local geometry optimization. In an AIRSS search, one can generate initial structures with a specified number of randomly chosen symmetry operations. The effect of symmetry will be discussed further in Sec. III A 1.

As a random sampling approach, there are no rigorous convergence criteria to tell when to stop an AIRSS search. However, it is clear when the search should not be stopped. Searching should continue at least until the following requirements are fulfilled: (1) known marker structures, if any, are identified; (2) a number of structures with low energies are encountered; and (3) these structures are encountered multiple times.

As part of an AIRSS search, the geometry optimization and the energy evaluations are carried out based on first principles DFT calculations. In this work, the plane-wave DFT code CASTEP was used for each geometry optimization.^36,37 The Perdew–Burke–Ernzerho (PBE) exchange–correlation functionals and the on-the-fly generated ultrasoft pseudopotentials (QC5) were adopted.^37,38 A relatively low plane-wave cutoff energy of 340 eV was used. The reciprocal space was sampled on Monkhorst–Pack grids with a spacing of 0.07 Å⁻¹. The effect of spin-polarization and the Hubbard U correction on the sampling of the PES will be discussed in Sec. III A 1.

For further structural optimization and property calculations, spin-polarized DFT calculations were performed using PBE exchange–correlation functionals with Hubbard U corrections. For this, we used the VASP code to make use of the calibrated U parameters for phase stability and voltage prediction.^39–43 A plane-wave cutoff energy of 520 eV was used, alongside Monkhorst–Pack grids with a spacing of 0.04 Å⁻¹. Hubbard U corrections of 5.3, 6.2, 3.25, 3.9, and 3.32 eV were applied to the d-channel of Fe, Ni, V, Mn, and Co, respectively.^41,43 The bandgap values were estimated using the HSE06 hybrid functional based on the PBE + U relaxed structures.⁴⁴ Ferromagnetic spin arrangements are assumed in spin-polarized DFT calculations.

The average cathode voltage V_a is calculated using the following equation:

where $ECat−Lix$ and E_Li are the enthalpy per atom of the cathode with a Li concentration of x and the enthalpy per atom of the lithium metal, respectively. e is the charge of an electron.

It is worth noting that the voltage of a cathode is relatively difficult to accurately determine from simple DFT calculations. One source of inaccuracy is the inexact exchange–correlation functionals employed in current DFT schemes.⁴⁵ While a number of studies have been carried out showing that high-level hybrid functionals give good approximations to the voltages, in this work, we chose to use the PBE functional with a Hubbard U correction, which is capable of predicting the voltages of cathodes with satisfactory accuracy and is sufficient for screening purposes.⁴⁶ Another source of inaccuracy in voltage prediction comes from an overestimated structural stability in DFT calculations. While in reality, the cathode may go through phase transition/structural degradation upon delithiation, DFT calculations are usually performed on small cells, and the relaxation scheme is not able to capture these transitions effectively. Therefore, the average voltage calculated by relaxing the delithiated structures is an estimate of the upper bound. Similar to the voltage, the capacity is difficult to be predicted accurately as well. In principle, the capacity is determined by the thresholds of Li concentrations below/above which the structure starts to collapse or degrade. However, large scale structural degradation is difficult to capture using DFT-based relaxation schemes as mentioned earlier. In this work, we estimated the upper bound of capacity by looking at the stability of the delithiated structures by structural relaxation. When over-delithiated, the C₂O₄ anions decompose into linear CO₂ molecules. Furthermore, the gravimetric energy density of the cathode is estimated by multiplying the average voltage and the capacity followed by normalization by its weight.

The rate capability of cathodes can be estimated by calculating the Li diffusion barrier and the bandgap, which give information about the ion and the electron transport, respectively. The Li diffusion barrier was estimated using climbing image nudged elastic band calculations,⁴⁷ and the bandgaps were calculated using the HSE06 functional as described previously.⁴⁴ 

To present the results of our searches, we generated structure maps depicting the relative positions of each structure on the underlying PES. To do this, we first generated the Smooth Overlap of Atomic Position (SOAP)^48,49 description of each structure in a given dataset using the ASAP code,⁵⁰ characterizing each structure as a vector in a high-dimensional space. We used the universal SOAP parameters in ASAP, which are constructed based on the elements present in the structures, following the heuristics of Cheng et al.⁵⁰ Note that we included crossover terms between elements in these descriptors, and the global description of each structure is yielded by averaging over the local environments on an element-wise basis. The descriptor vectors were normalized so that each variable had a variance of unity across a given dataset.

We then applied the dimensionality reduction method Stochastic Hyperspace Embedding And Projection (SHEAP) to the high-dimensional data to produce two-dimensional representations.^51,52 SHEAP constructs a weighted graph describing the similarity relationships in the source data. Following a test to combine equivalent structures, the data are cast into a low-dimensional space (2D in this case) using random projection. Another graph of weights is then constructed for the nodes defined by these mapped points. The layout of the map is then optimized by minimizing a cost-function that penalizes any mismatch between the weights.

In these maps, individual structures are represented by circles colored according to enthalpy, with areas proportional to the number of occurrences in the search. Circles are prevented from overlapping each other by a soft-core repulsion that is turned on only once the nearly optimal layout has been reached. The axes do not have any predetermined physical interpretation, but the relative positioning and clustering across the two-dimensional space reveals the structural similarity.

We start by performing an AIRSS search on LiCoO₂, the most commonly used cathode material in LIBs. For this system, we conducted the search using the minimum separation matrix and volume per formula measured from the known $R3¯m$ high-temperature (HT) phase.^53,54 The values are listed in Table I. The number of formula units was chosen to take values of 1–5. The results of this search are summarized in Fig. 2(a). In total, 1200 random structures were generated and optimized. Of the relaxed structures, the trigonal-layered HT phase had the highest encounter rate at 47%. The cubic-lithiated-spinel low-temperature (LT) phase, which has a comparable energy, was also found multiple times (∼2%).⁵⁵ Other low energy structures found include the P2 and O2 phases of LiCoO₂, which result from a slight variation of the CoO₆ layer stacking compared to the O3-type stacking of HT-LiCoO₂.^56,57 All of these low energy structures are characterized by similar local coordination—cobalt ions are octahedral, while lithium ions are either octahedral or prismatic. This search also turned up many polymorphs with relatively high energies. These include several approximants of disordered rocksalts, in which cobalt and lithium ions are distributed in a non-layered manner, within a distorted MO₆ (M = Li and Co) framework. In these structures, the large lattice distortion gives rise to the higher energies. We expect these disordered rocksalts to be difficult to synthesize.

The above-mentioned search was carried out based on prior knowledge of LiCoO₂. However, in exploratory searches, such prior knowledge may not exist. Now, we demonstrate that constraining the search space using the minimum separation matrix, cell volumes, and symmetry is critical to achieving efficient sampling. This is followed by a discussion on how one can use a preliminary AIRSS search to determine a suitable list of minimum separation values and cell volumes.

Figures 2(b) and 3(c) summarize the search results obtained from the use of a matrix of minimum separations drawn uniformly from 1–3 Å and of a single fixed minimum separation for all species pairs of 1.5 Å, respectively. In comparison to Fig. 2(a), these searches resulted in many more high energy configurations with a significantly lower encounter rate for the lowest energy structures. In the case of a single fixed minimum separation, we can easily understand this by noting that the chemical environments of the cation and anion are not distinguished by this constraint. As a result, Li–Li and O–O close contacts are generated resulting in unphysical and high energy structures and chemistries.

As illustrated by the results of Fig. 2, an effective way to select suitable minimum separation values is to measure them from known structures with similar chemical compositions. However, in some cases, the detailed coordination environment is not known as a prior, and the minimum separations need to be selected using a first principles approach. A way around this is through a preliminary AIRSS search. Figure 3 shows the distribution of minimum separations arising in relaxed structures obtained from a search on LiCoO₂ using a random minimum separation matrix whose elements are drawn uniformly from 1–3 Å. No symmetry constraints were applied. We observed, for each pair of species, rapid convergence in the lower bound of the interatomic distances with respect to increasing the number of formula units, with the converged values corresponding to the minimum bond lengths of the system. Thus, for a “safe” search, one can use a relatively small cell to acquire these minimum separations, ruling out certain close contacts between pairs of species. However, these lower bounds do not consider the detailed local coordination of the system. For example, as shown in the inset of Fig. 3, the lower bound of Co–Co distance is ∼2.1 Å, which corresponds to Co pairs in an O-connected CoO₄–CoO₄ configuration. Compared with the CoO₆ octahedra, such local coordination is energetically less favorable. Therefore, for an “efficient” search, it is better to choose the minimum separations from the structure with the lowest energy, as indicated by the dashed lines in Fig. 3. This will bias the search toward having local coordinations that are more likely to occur in low energy structures.

A search can be further speed up through the application of symmetry constraints. As illustrated in Figs. 2(a) and 2(d), for a well-chosen minimum separation matrix, the encounter rates for low energy structures are comparably high with and without the use of symmetry. However, the overall searching time is significantly longer (∼10 times) in the absence of symmetry due to slower DFT electronic steps during the geometry optimizations.

It is worth noting that the majority of the time taken for an AIRSS search is spent on the first principles geometry optimizations. Thus, speeding up the DFT calculations is highly beneficial. One way we achieved this was by adopting ultrasoft pseudopotentials (QC5) specially made for high-throughput calculations, which allow the use of low energy cutoff. Apart from that, the inclusion of spin polarization and Hubbard U parameters also affects the speed of calculation. Therefore, when necessary and tested, another way to accelerate the process is to perform non-spin polarized calculations during the search followed by a subsequent round of high accuracy calculations to refine the low energy structures. For LiCoO₂, we compare the structural density of states obtained without spin, with spin, and with the Hubbard U correction. Here, the structural density of states is described as the proportion of structures at each energy. It evaluates the frequency of finding a local energy minimum with a specific energy on the potential energy surface. We found that the structural density of states obtained using different computational parameters is similar, and the low energy structures are all found with comparable encounter rates, as shown in Fig. 4(a). Figure 4(b) further shows the correspondence of energies obtained with non-spin polarized PBE and with spin-polarized PBE + U for the structures of LiCoO₂. We observe a shift of stability order with an energy variation of ∼150 meV/atom.^32,33 For an AIRSS search carried out using non-spin polarized PBE, any structure with energy below this value would need to be recalculated to obtain the correct ordering. The magnitude of this threshold is of course system dependent. In extreme cases, the choice of DFT parameters may significantly alter the PES, possibly resulting in missing certain local minima. Therefore, we must note that computational parameters (e.g., exchange–correlation functionals, spin-polarization, and correction schemes) must be chosen with caution. Preliminary and post-mortem searches need to be carried out to check whether the results are parameter dependent. If they are, the best practice is to report a confidence interval instead of a single prediction.

In summary, our workflow for efficient AIRSS is as follows: First, obtain the species-pairwise minimum separations and the cell volume parameters, which constrain the search space, either by considering prior knowledge of the chemical system or by carrying out a preliminary search using random minimum separations; then, generate and relax random sensible structures with different numbers of formula units. Particularly in the case of many formula units, symmetry may need to be taken into consideration to accelerate the search. Additional re-optimization of low energy structures may need to be carried out if the search used relatively strong assumptions for the DFT calculations.

LiFePO₄ is a representative polyanion cathode material that is widely used in electric vehicles due to its low raw material cost.⁵⁸ LiFePO₄ is structurally more complicated than LiCoO₂ and therefore is a more challenging task for crystal structure prediction. In our application of AIRSS to this system, we generated random structures fulfilling the following requirements: (1) the structure has 2 to 4 symmetry operations; (2) the P and O species are generated as proper PO₄ units; (3) the interatomic distances are larger than the predefined species-pairwise minimum separation values, which are obtained from the known olivine phase (the P–O distances are specially treated so that only those measured from inter-PO₄ units are constrained to fulfill this requirement); (4) the volume per formula unit is within 15% as compared with that of the known olivine phase. The search results obtained from relaxing these random sensible structures are summarized in Fig. 5. Out of the 684 unique structures found, the experimentally observed olivine phase (Pnma) turned up four times, giving an encounter rate of ∼0.58%. Also found was a Cmcm structure, slightly higher in energy than the olivine phase, which has previously been reported under high pressure conditions.⁵⁹ Interestingly, many other phases with similar energies to this were located. Some of these are likely to be synthesizable under fine-tuned conditions and thus may be worth further computational or experimental study.

In our study of this system, we observed that the use of pre-defined chemical subunits (PO₄), as opposed to individual atoms (P and O), was effective in steering the search toward lower energy structures. This is illustrated in the inset of Fig. 5 through the structural density of states obtained with and without this constraint. Caution is needed with such constraints, however, since they restrict the sampling space. They should be used only when the chemistry is well understood.⁶⁰ In the case of LiFePO₄, structures containing other polyphosphate anions are not energetically competitive and so can be disregarded, allowing us to accelerate the search by imposing P and O to be present as tetrahedral units of PO₄.

The cases of LiCoO₂ and LiFePO₄ demonstrated the power of AIRSS for compounds of a given single composition. However, in exploratory searches, one does not know in advance which composition will give the best performance or the desired properties. Therefore, multi-compositional searches followed by post-search screening are necessary. More importantly, the determination of the stability of an unknown phase not only involves the calculation of its formation energy but also whether it is susceptible to decomposition into other phases. Therefore, it is necessary to predict the entire compositional phase diagram. For cathode materials, this is especially important when the chemistry of the system is less understood. The construction of the phase diagram is critical in determining the reaction pathways during lithiation since not all cathodes are of the intercalation type and topotactic phase transformation may happen. For example, transition metal fluorides, for example, FeF₂, FeF₃, and CuF_2, are typical cathode materials that go through a conversion reaction during discharge.⁶¹ Here, we consider the system of Li–Cu–F, demonstrating how a ternary phase diagram can be generated by conducting an extensive multi-composition search in the chemical space.⁶² In particular, we searched for structures with compositions with the formula Li_xCu_yF_z, where x, y, and z fulfill the following requirements:

and

where $Z+$ is the set of positive integers. This led to 163 compositions. For each composition, we generated at least 50 random structures for further first principles optimization. In total, 15 031 structures were generated and optimized. This number was found just enough to cover the stable phases in this study. In exploratory searches where the chemistries are unknown, many more structures may need to be considered. The number of formula units was restricted to being 1 and 2 in the current search. The minimum separation matrix and the cell volumes were determined from a preliminary search that used random species-pairwise minimum separations of 1–3 Å with ∼50 structures being sampled on each composition. The details are discussed in Sec. III A 1. All the structures generated were constrained to have 2 to 4 symmetry operations.

As shown in Fig. 6(a), multiple stable phases were obtained from the search. For the binary compounds, LiF, CuF₂ and Cu₃Li were found to be stable. Such a result helps determine the terminal redox couple of CuF₂ as a cathode, i.e., the final thermodynamically stable reaction product should be Cu and LiF. The average voltage is calculated to be 3.64 V. For ternary compounds, multiple compounds were found to shape the convex hull, including Li₂CuF₆, LiCu₂F₆, and LiCuF_4, in close agreement with existing databases.¹⁶ Since none of these stable phases are located within the CuF₂–LiF–Cu triangle, CuF₂ is expected to adopt a single-step conversion reaction to form LiF and Cu upon lithiation. Interestingly, we also found a number of phases that are less stable (above the convex hull) but are potentially interesting based on kinetics. For example, the phase stability along the Li–CuF₂ tie line is shown in Fig. 6(b). The hypothetical Li-intercalated compound LiCuF₂ is high in energy (220 meV/atom above hull). Such an instability explains why CuF₂ goes through a conversion reaction, rather than intercalation, when lithiated. We also found phases containing F⁻ anion lattices that are topotactical with those in LiF, as shown in Fig. 6(c). These have recently been proposed as possible buffer phases between the reactant CuF₂ and the reaction product LiF.⁶³ Cu (I) containing fluorides have also been proposed as potential intermediate phases, especially during the charging process. Xiao and co-workers proposed Cu (I) containing phases with tetrahedrally coordinated Cu following electron diffraction analysis.⁶² Our search reveals possible phases that are consistent with this. As shown in Fig. 6(e), the Li₃CuF₄ phase has an energy above the hull of 60 meV/atom. Despite its relatively high energy, the formation of this at interfaces and under electrochemical conditions, is plausible.

To conclude, a multi-composition search can efficiently generate the phase diagram for prediction of phase stabilities and for analysis of the reaction mechanisms of a cathode.

Transition metal oxalates are a family of polyanionic compounds that possess a number of advantageous characteristics as cathodes for batteries but have been overlooked as compared to other polyanionic materials such as the phosphates PO₄³⁻, sulfates SO₄²⁻, and silicates SiO₄⁴⁻.^64,65 Interestingly, the C₂O₄²⁻ anion has a similar polarizability to PO₄³⁻. As a result, it should be able to deliver reasonable redox potentials during lithiation, such as LiFePO₄. Furthermore, oxalates can be synthesized under mild conditions via solution-based methods. In fact, they are known to nucleate into a chain-like, layered, and fully connected frameworks, thanks to complex ligand coordination chemistry. This offers exciting opportunities for oxalates to serve as cathode materials by tuning the structural features through metastability. Despite these advantages, oxalates have been overlooked as cathode materials, in part, due to the large molecular weight of the anion, which leads to relatively low energy density.⁶⁵ Recently, with the discovery of the reversible anion redox as an additional source of capacity, the oxalates are gaining attention from battery researchers. Jiao et al. synthesized a lithium-containing iron oxalate with composition Li₂Fe(C₂O₄)₂.⁶⁶ In this material, Fe has an oxidation state of +2. Interestingly, they found that more than 1 Li/f.u. can be taken out of the structure reversibly. Using spectroscopic evidence, they found that this over-delithiation is enabled by the reversible electron stripping from the C₂O₄²⁻ anion. As compared to the oxides, especially Li-rich ones, the strong covalent O–C bonds in C₂O₄²⁻ prohibit the O₂²⁻-dimer formation and the O₂ evolution during anion oxidation.

Given the possible polymorphism and interesting chemistry in Li-stuffed transition metal oxalates, we carried out exploratory AIRSS searches on a range of systems of this nature. In particular, we considered Li₂TM(C₂O₄)₂, where TM = Fe, Co, Ni, V, and Mn. The aim was to find low energy phases with the potential to be synthesized and display reasonable energy density and rate capability as cathodes for LIBs.

Figure 7(a) summarizes our AIRSS searches for Li₂Fe(C₂O₄)₂ in the form of a SHEAP map. During the search, C and O were treated as C₂O₄ units when generating the initial random structures, and the number of formula units was restricted to being 1–4. As illustrated in the SHEAP map, among the large number of local energy minima found, only a fraction are found to be low in energy, emphasizing the complexity of this PES. Representative low energy structures are labeled in Fig. 7(a). The structures with the highest encounter rates are layered, in which the Li⁺ cations are distributed within the layers. These structures resemble the high-temperature phase of LiCoO₂. A number of “checkerboard” structures, composed of one-dimensional structural units of octahedral-coordinated Fe and planar square-coordinated Li, were also encountered frequently.

Table II lists the energetics of the low energy phases of Li₂TM(C₂O₄)_2, whose structures are given in Figs. 7(d)–7(i) using Fe compounds as an example. Interestingly, for Fe, Co, and Ni, we found that at least one of the layered polymorphs has an energy lower than the experimentally reported structure at the PBE level. For Fe, this polymorph is characterized by its well-connected Li-network and is expected to conduct lithium more efficiently than the experimental phase.⁶⁷ This is confirmed through calculation of the diffusion barrier of Li in the two structures. In the layered structure, the barrier for Li diffusion is ∼320 meV, close to some superionic conductors, whereas for the experimentally reported structure, the value is ∼620 meV.^66,68 This indicates that the layered phase has the potential to display high-rate capability as a cathode.

In our search, as well as this layered ground-state structure, a number of metastable structures were found that may be synthesizable through complex nucleation kinetics in solution. Indeed, we calculated that the experimentally reported Li₂Fe(C₂O₄)₂ phase is also intrinsically unstable, with an energy above the hull of 58 meV/atom. Here, we used such a value as a reference and increased it slightly (by 20 meV/atom) as the threshold for selecting Li₂TM(C₂O₄)₂ structures on which to conduct further property calculations. These are the structures before the dashed line in Fig. 7(c). The average discharge voltages of structures below this threshold are listed in Table II. For Fe, all of them are over 4V vs Li/Li⁺, as expected, since the anion redox is involved. We estimated upper bounds for capacities of these materials by checking the structural integrity of the delithiated phases at different states of charge. We found that when the structure is over-delithiated, the C₂O₄²⁻ anion decomposes into two CO₂ molecules. Using whether such decomposition happens during structural relaxation as a criterion, the estimated capacity upper bound is listed in Table II. Interestingly, among the TMs we considered, V-oxalates seemed to give the highest capacity without structural collapse, whereas Ni-oxalates easily decomposed when delithiated. The stability of V-oxalates may be due to the multi-valence nature of V. Furthermore, by multiplying the capacity and the voltage, we computed an estimate for the energy density of oxalates. A number of polymorphs (Fe-Layered #1, Mn-Layered #2, Mn-Mixed, Co-Layered #2, and all V-polymorphs) appear to have the potential to reach extraordinary gravimetric energy densities (>900 Wh kg⁻¹), close to or even higher than that of LiNi_0.8Co_0.1Mn_0.1O₂. It is worth noting that, in these predictions, temperature effects are not considered. The entropy term arising from partially occupied Li sites/different Li arrangements at finite temperatures is not incorporated. While this will not significantly alter the energy densities predicted here and is sufficient for screening purposes, identifying disorder and computing the corresponding free energy is critical in analyzing the phase separation behavior of cathodes during discharge. Further methodological development is required in this context, such as the incorporation of lattice-based sampling methods, for example, cluster expansions.

We have explored ab initio random structure searching as a tool for battery cathode discovery. We showed, through case studies of known cathode materials of LiCoO₂, LiFePO₄, and Li_xCu_yF_z, that efficient prediction of the low energy structures can be achieved. We demonstrated that using species dependent minimum interatomic separations and cell volumes to steer the search is critical to its efficiency. Further acceleration was achieved by imposing symmetry constraints and by using pre-defined chemical subunits for structure generation. Based on these principles, we have carried out explorative searches for Li₂TM(C₂O₄)₂ oxalates. A number of low energy polymorphs with layered structures were found. Post-search screening and property calculations identified promising cathode materials characterized by high voltages (>4V) and low Li diffusion barriers (∼300 meV), with the potential to deliver not only high energy densities (up to ∼900Wh kg⁻¹) but also high-rate capability. This work serves as a basis for AIRSS-based cathode searches, providing a detailed framework for an efficient workflow.

Z.L. and B.Z. contributed equally to this work.

This work was supported by the Faraday Institution (Grant No. FIRG017) and used the MICHAEL computing facilities. C.J.P. acknowledges support from the EPSRC through the UKCP consortium (Grant No. EP/P022596/1). D.O.S. acknowledges support from the European Research Council, ERC (Grant No. 758345). Through the membership of the UK’s HEC Materials Chemistry Consortium, which is funded by the EPSRC (Grant Nos. EP/L000202 and EP/R029431), this work used the ARCHER UK National Supercomputing Service (www.archer.ac.uk) and the UK Materials and Molecular Modeling (MMM) Hub (Thomas Grant No. EP/P020194 & Young Grant No. EP/T022213). B.W.B.S. acknowledges EPSRC CDT in Computational Methods for Materials Science for funding under Grant No. EP/L015552/1.

The data that support the findings of this study are available from the corresponding author upon reasonable request.