Low-cost descriptors of electrostatic and electronic contributions to anion redox activity in batteries

Conventional battery cathodes are limited by the redox capacity of the transition metal components. For example, the delithiation of LiCoO2 involves the formal oxidation from Co(III) to Co(IV). Enhanced capacities can be achieved if the anion also contributes to reversible oxidation. The origins of redox activity in crystals are difficult to quantify from experimental measurements or first-principles materials modelling. We present practical procedures to describe the electrostatic (Madelung potential) and electronic (integrated density of states) contributions, which are applied to the LiMO2 and Li2MO3 (M = Ti, V, Cr, Mn, Fe, Co, Ni, Cu, Zr, Nb, Mo, Ru, Rh, Pd, Ag, Hf, Ta, W, Re, Os, Ir, Pt, Au) model systems. We discuss how such descriptors could be integrated in a materials design workflow.

The large difference in electronegativity between metals and oxygen mean that metal oxides are heteropolar compounds [1,2]. Oxidation of metal oxides occur between two extremes. In the first, common for transition metal oxides, an increase in the formal oxidation state of the metal occurs ( ), resulting in a so-called oxygen hole [3]. Due to the electron configuration of the oxide ion going from closed-shell 2p 6 to open-shell 2p 5 , evidence for oxygen hole formation has commonly being reported from electron paramagnetic resonance (EPR) spectroscopy. As early as 1955, O'Brien and Bleaney used EPR to establish that the colour centres in smokey quartz (SiO 2 ) are due to oxygen redox induced by Al impurities [4]. EPR signals are presented as clear evidence for the hopping of holes between oxygen sites in SiO 2 in the later work of Schnadt and Schneider [5].
The dominant oxidation mechanism in conventional intercalation battery cathodes is cationic, while recently there has been a strong focus on the development of so-called anion redox cathodes. In this context, anion redox refers to oxidation processes in which M and O are both active [6], which could enhance overall battery capacity. Detailed reviews on this topic can be found at [7] and [8].
The driving forces dictating to what extent M and O can be oxidised in the same compound are challenging to probe by experiment and simulation. Recently, there have been substantial efforts made to understand how anion redox might be harnessed in a small set of emerging cathode materials [9][10][11][12][13][14][15][16].
Meanwhile, developments in predictive, high-throughput design in other areas of energy materials research such as photovoltaics [17][18][19], photocatalysis [20,21] and thermoelectrics [22][23][24] have enabled broader searches of the chemical landscape and have in some cases led to the discovery of promising candidate compounds. Such virtual screening studies rely on the use of simple descriptors to predict properties of candidate materials at low computational cost. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
Here we present practical procedures for probing the contributions to redox phenomena from electrostatic and electronic perspectives. Our workflow relies on open-source materials informatics tools and the output of routine first-principles calculations. We provide worked examples to demonstrate our approach, which is applied to LiMO 2 and Li 2 MO 3 model systems (M=Ti, V, Cr, Mn, Fe, Co, Ni, Cu, Zr, Nb, Mo, Ru, Rh, Pd, Ag, Hf, Ta, W, Re, Os, Ir, Pt, Au), that have been structurally optimized using density functional theory (DFT).

Electrostatic descriptors
Oxide (O 2− ) ions are unstable in vacuum, but become stable in a crystal due to the confining electrostatic potential [25]. This makes the positioning of the O 2p valence band in a metal oxide highly sensitive to the crystal environment [26].
The electrostatic (Madelung) energy of a crystal under the ionic point-charge approximation can be expressed as the sum of all pairwise interactions, where z i and z j are the formal charges of the ith and jth ions and r ij is the interionic distance. Accordingly, the potential experienced at site i (the site potential) is . A larger potential is produced for more highly charged metals at shorter distances.
A smaller average Madelung site potential might be expected from a higher Li + :M n+ ratio, owing to the overall decrease in positive charge. For stoichiometric systems, however, the value of n increases with this ratio in order to satisfy the valence of all species, and this is the dominant effect as seen in Figure 1 when comparing LiMO 2 (M 3+ ) to Li 2 MO 3 (M 4+ ) for a given transition metal.
There is a range of >3.5 V for the average Madelung potential of oxygen even within a given structure prototype, (LiAuO 2 vs LiCoO 2 , Figure 1(a)). Given that the formal charge of the transition metals remains constant within a series, the variation in V i is solely to differences in interionic distances in the relaxed crystal structures. It should be noted that in addition to the compounds containing group 10 metals (Ni, Pd, Pt), which relax to a different structure type with square planar geometry, Nb, Ta and Os also distort the parent structure shown in Figure 1.
If Li/M site disorder is introduced-based on 128 atom 2×2×2 supercells of LiMO 2 -the impact on the local site potentials is large. This effect is apparent in Figure 2, which shows the oxide potential range in LiMO 2 for the first five transition metals that we consider. The site potential decreases rapidly as the Li + :M 3+ ratio increases in a similar way for all five compounds. This trend recovers the analysis of Seo et al, which showed that a Li-rich coordination is required for the extraction of oxygen electrons in Li-ion cathodes [11].  [27]. A smaller value of DE h indicates more favourable hole formation on oxygen. Figure 3 shows how this descriptor largely follows the same trend as the negative of I n M , but is modulated by the structural contributions via the Madelung potential. Note that this trend is most obviously observed for the first row transition metals (Figure 3) for which a complete set of reliable ionisation energies is available.
In summary, the Madelung site potential is a physically meaningful metric that is simple to calculate. While there are few theoretical or experimental results with which the observations made on our model systems can be directly compared, we have shown how the O site potential can be highly sensitive to the surroundings in typical Li-ion cathode compounds, that the observed trends are consistent with recent theoretical work in the field of anion redox [11], and how this can be extended to estimate hole localisation energy.

Electronic descriptors
A quantitative description of cathode redox processes require explicit consideration of delithiation. In the spirit of finding computationally cheap descriptors, however, we can examine the lithiated material [6,7].
The electronic density of states (DOS) represents the number of states available as a function of energy, and is a standard output of all electronic structure codes. We can project this onto atomic sites to obtain a species and orbital decomposed DOS. By integrating the oxidation window in the upper valence band, it is possible to extract  the ratio of metal d-orbital (M d) to oxygen p-orbital (O p) contributions. Here we choose an oxidation window of 1 eV below the top of the valence band, but this could be easily adapted to more accurate distributions. Figure 4 shows the extreme cases of MgO and Cu 2 O, in which Cu + has a d 10 configuration. DOS analysis (Notebook 2) shows that the upper valence band is comprised of 100% O p and 0% Mg s for MgO, and 87% Cu d and 13% O p for Cu 2 O. These fractions are expressed in terms of the orbitals of interest (e.g. Cu d + O p=100% for Cu 2 O) to avoid influence on the stoichiometry, which will change upon cycling.
The DOS fractions of M-d and O-p are compared for LiRuO 2 , Li 2 RuO 3 , LiWO 2 and Li 2 WO 3 in Figure 5. The O-p contribution to the upper valence band is higher for Li 2 MO 3 compared to LiMO 2 . It should also be noted that the ground state magnetic ordering of Li 2 RuO 3 is ferromagnetic (Figure 5), which is necessary to obtain a reliable electronic structure. An efficient way to automate the search for magnetic ground states is described by Horton et alin [29].
This analysis allows for an immediate insight into which states are available in the upper valence band to contribute to redox processes. Our approach exemplifies how tracking individual orbital contributions to the DOS, which has constituted important evidence for the presence of absence of anionic redox in previous studies [12,30], can be quantified.

Workflow integration
The electrostatic descriptors are low-cost and only require structural information. Calculation of Madelung site potentials can be carried out in a few seconds on a desktop computer using the Ewald method, which is implemented in Pymatgen [31] (see Notebook 1). While Madelung energies have previously been used as descriptors of battery voltage [32], individual site potentials have been applied to understanding band  alignments [33,34]. The approaches could be improved by incorporating dielectric polarisation of the host crystal.
The electronic descriptors require only the DOS, and this is produced as standard output for DFT calculations. The post-processing is performed in less than a second on a desktop computer; a negligible amount of time compared to the original DFT calculations (see Notebook 2).
Overall, both descriptors require common representations of crystalline compounds (i.e. crystal structure and electronic DOS) and are therefore suitable for large-scale studies of lithiation, sodiation in Na-ion battery cathodes, and other chemistries. They could both be integrated, for example, into Monte Carlo techniques as part of multi-scale simulations ( Figure 6).

Technical details
All the code and data necessary to reproduce the results discussed here can be found the accompanying data repository at https://doi.org/10.5281/zenodo.3754712. We rely primarily on the Pymatgen Python library [31], which is used to process crystal structures optimized by DFT (Notebook 1) and electronic structures solved by hybrid DFT (Notebook 2).
All first-principles calculation were automated using the Python packages Atomate [36] and Fireworks [37]. We used Kohn-Sham DFT with a projector-augmented plane wave basis [38] as implemented in the Vienna Ab-initio Simulation Package (VASP) [39,40]. We use the PBEsol exchange-correlation functional [41] and a k-point grid is generated for each calculation with a density of 120 Å 3 in the reciprocal lattice. The kinetic-energy cut-off is set at 600 eV and the forces on each atom minimised to below 0.005 eVÅ −1 .
Semi-local exchange-correlation treatments such as the PBEsol functional provide an accurate description of crystal structures but tend to underestimate the electronic bandgaps of semiconductors. To overcome this issue, more accurate calculations are performed using a hybrid non-local functional HSE06 [42], which includes 25% screened Hartree-Fock exact exchange, which are used for the electronic structure analysis. Γ-centred homogeneous k-point grids are used with a density of 120 Å 3 in the reciprocal lattice and the kinetic energy cutoff is set at 550 eV.