DFT calculations have become an important part of materials research to discover and explain the causes of experimentally observed phenomena at the atomic scale. They provide insights into the physics and chemistry of materials which aid in further optimisation of materials for a specific application. DFT primarily provides a means for calculating the total energy and electron charge distribution of any configuration of atoms.

The atomic-scale processes in materials are described by the quantum mechanical time-independent Schrödinger equation, Equation (i):

(i)

in which the wavefunction for the set of electrons and nuclei is denoted by Ψ({R_j},{r_i}) where R_j are the positions of the nuclei, r_i are the positions of the electrons and Ĥ is the Hamiltonian of the system. The energy E obtained from this equation represents a specific energy level for the system. In general, the ground-state energy of the system, E₀, is the quantity of interest. The Hamiltonian for this time-independent equation is Equation (ii):

(ii)

The first two terms in Ĥ are the kinetic energy operators of the nuclei and electrons, and the third is the potential energy. Nuclei and electrons interact via the Coulomb interaction. Unfortunately, the conventional Schrödinger equation is too complicated to solve beyond just a handful of particles. Therefore, approximations are required in order to solve this equation and obtain the ground-state energy of the system of interacting electrons and nuclei. Since electrons move on very fast timescales compared to nuclear motion, the nuclei can be treated as fixed in space while the electronic-ground state is computed. This is the Born-Oppenheimer approximation, which results in a Schrödinger equation for the electrons, in which the nuclear positions and charges enter as parameters only. The underpinning principle of DFT, the Hohenberg-Kohn theorem (1), builds from this approximation, providing a theoretical basis for working not with the wavefunction, but with the much simpler ground-state electron density, n(r).

Figure 1 shows an example of the calculated ground-state electron density of the atoms in a silicon crystal structure, represented by the smooth surface surrounding the atoms. The total energy of a system of electrons and fixed nuclei is a function of all possible electron density functions. Using the Kohn-Sham ansatz, finding the ground-state electronic density is made computationally feasible by expressing it in terms of auxiliary wavefunctions which describe a fictitious non-interacting system of the same density (2). The full expression for the ground-state energy E_KS may then be written as Equation (iii):

(iii)

Fig. 1.

where the first term, T[n], is the kinetic energy associated with the non-interacting Kohn-Sham particles; the second term, E_NN, is the nuclear-nuclear interaction; and the third term, V_ext, is the external potential of ion cores in which the electrons move. The fourth and fifth terms represent electron-electron interaction energies. The fourth term is the exact classical electrostatic energy; the interaction energy of an electron with the mean field of all electrons. The fifth term is the exchange-correlation energy, which attempts to account for all interactions not accounted for within the first four terms. By dividing up the energy in this way, while the exact exchange-correlation functional remains unknown, it may be approximated in various tractable ways.

The simplest approximation to E_XC is the local density approximation (LDA), where the exchange-correlation energy per particle is taken to be equal to that of a uniform electron gas of the same electron density, at each point in space. Generalised gradient approximation (GGA) functionals improve on the LDA by taking into account both the electron density and the gradient of that density, resulting in a more accurate description of exchange and correlation (3). These functionals have limitations; most seriously, both electron localisation and electronic band gaps are underestimated. So‐called ‘hybrid functionals’ have aimed at semi-empirically correcting the electronic band gap (4) and developing functionals beyond the LDA and GGA is the focus of much of the theoretical work in the field of DFT today, where the ultimate goal is to find an exchange-correlation functional which accurately describes all possible systems (5).

Within this framework, total energies, forces, equilibrium geometries, elastic behaviour and many other properties of interest can be readily and accurately predicted. However, to predict a material’s properties using DFT, it is necessary to know how its atoms are arranged. Thus, in the following section, we describe the method of CSP, which uses DFT to generate structures of novel materials.

There are multiple materials databases. Some contain only the experimental crystal structures and other relevant properties of known materials, while others contain the computed properties of both known and hypothetical materials. These can be leveraged to perform CSP. For example, known crystal structure prototypes can be decorated with any set of atomic species, resulting in new hypothetical materials. The stability and synthesisability of these new materials can then be assessed using DFT calculations and by comparing against thermochemical data in the database. Three of the major exhaustive databases of DFT calculations, the Open Quantum Materials Database (OQMD) (6), the Automatic Flow (AFLOW) framework for materials discovery (7) and the Materials Project (8) have been used to predict new materials and screen for desired properties using a combination of high-throughput ab initio calculations and, increasingly, statistical and machine learning approaches. In addition, experimentally identified structures are found in the Inorganic Crystal Structure Database (ICSD) (9) and the Crystallography Open Database (COD) (10). These databases have been used as a starting point for many theoretical studies, leading to several new discoveries in the field of energy storage, including identifying SrFeO_3-δ as a material for carbon capture (11), verifying Li₃OCl as a solid electrolyte with high ion conductivity (12) and predicting LiMnBO₃ as a Li-ion battery cathode (13). While these databases are useful for comparisons of known structures and enable the discovery of materials that are based on known crystal structure prototypes, it is likely that new structures exist which cannot be classified as one of the currently known prototypes. Therefore, it is necessary to perform CSP in order to explore novel phases of materials.

The search for new thermodynamically stable materials (those favoured to form during synthesis, when kinetic factors are excluded) using CSP can take one of many approaches (14), but all involve a search for the lowest energy minimum in a high-dimensional configuration space. The configuration space for a periodic structure with N atoms per unit cell has dimension 3N+3, taking into consideration the rotational symmetries and unit-cell degrees of freedom, whilst the number of local minima in the space scales exponentially with N (15). Ideally, all low-lying minima would be sampled during CSP since metastable phases may be synthesised experimentally, or indeed be thermodynamically stable under different conditions; for example, graphite is the most stable allotrope of carbon under ambient conditions, but diamond can be easily synthesised under high pressure. Particularly popular approaches to CSP include the use of evolutionary algorithms to ‘breed’ new structures (15) and particle swarm optimisation (16–18).

AIRSS (19) is the focus of this review. Despite the potential for having a high computational cost, AIRSS remains an effective method for structure prediction which allows for a breadth of searching and has proven successful in a wide range of materials. Beyond the ease of its implementation, AIRSS has several advantages. Firstly, individual relaxations do not depend on one another, hence all trials can be run concurrently making the algorithm trivially parallelisable to the largest of supercomputers. Secondly, AIRSS allows for the easy application of chemically intuitive constraints which reduce the initial search space to the most experimentally relevant trial structures. This constraint greatly reduces the size of the search space and makes AIRSS applicable to a wide range of systems, including those at high pressure (20, 21). These chemical constraints include, for example: the phases of conversion and alloying anodes (22–25) were constrained by space group symmetries and atomic distances; high pressure phases of ice (26) were constrained to H₂O units; encapsulated nanowires (27) were constrained by rod group symmetries; metal-organic frameworks (28) were constrained to molecular building blocks; grain-boundary interfaces (29) and point-defects (30) had some atoms fixed to describe the lattice and systematically randomised other atoms to describe interface and defect structures.

AIRSS explores configuration space using random sampling as shown in Figure 2(a) and proceeds as follows.

Fig. 2.

To search for a new phase with chemical formula, A_x B_y , any number of atoms of element A and B are placed randomly (denoted ‘Randomise’ in Figure 2(a)) into a 3D simulation cell in the ratio x :y . The cell and atomic positions are allocated such that they obey a set of chosen symmetry operations (a space group in 3D). Further constraints, such as minimum separation between atoms and a feasible range for the atomic density of the unit cell, may be imposed. These constraints narrow the region of the configuration space of possible structures by avoiding regions that describe unrealistic arrangements of atoms.

The forces on the atoms and stresses on the cell are calculated with DFT and then minimised using the traditional optimisation algorithms (for example, conjugate gradients). This step is denoted ‘Relax’ in Figure 2(a). The energy of the system is used as a metric to gauge how stable the structure is.

Steps 1 and 2 are then repeated several thousand times in order to generate a representative set of structures in the A-B chemical space. The search is stopped once the lowest energy structures have been found multiple times. The set of lowest energy structures are the candidates for phases that are likely to form experimentally.

Using the DFT energies, one can construct a ‘convex hull’ of all the structures found by AIRSS, as shown in Figure 2(b). The structures A_x B_y which are likely to form, must both have a negative formation energy relative to elemental A and B and lie on the convex hull tie-line between A and B to avoid decomposition into other binary phases. This tie-line is shown by the black line connecting the lowest energy structures in Figure 2(b). This figure illustrates the process of constructing a convex hull using the optimised structures from AIRSS. Suppose at a given point during the AIRSS search, the only structures on the convex hull are AB, AB₂ and A₂B, connected by the dashed line in Figure 2(b). Subsequently, a novel phase, A₃B, is identified using AIRSS and is found to lie below the existing tie-line. In this case, CSP has identified a new ground-state structure which suggests an additional phase, A₃B is likely to exist within the A-B phase diagram. Therefore, as shown in Figure 2(b), the hull is reconstructed to include the phase A₃B, rendering A₂B unstable, given that it is now no longer on the convex hull. Although this example is given for two dimensions (i.e. a binary system containing elements A and B) the convex hull construction is generalisable to N dimensions, in which the tie-lines between the lowest energy structures are computed in a similar manner.

In this way, AIRSS enables the prediction of new thermodynamically stable and metastable compounds in a given phase diagram and the convex hull construction provides a guide to their stability compared to previously known phases, without performing exhaustive chemical synthesis. Synthesis experiments can then be targeted at the most promising compositions and characterisation experiments can be guided by the predicted model structures.

The electrochemical voltage profile is the voltage signal of the electrode measured (vs. a reference, usually Li⁺/Li) as a function of the number of ions (i.e. charge) stored in the electrode. The phase transitions, which occur within the electrode during cycling provide the characteristic shape of the voltage profile; two-phase regions show a constant voltage, while solid-solution regions show a sloping voltage. The voltage drop between two phases is proportional to the difference in their free energies and thus these voltage drops can be computed directly from the free energies of the phases which lie on the convex hull tie-line. The voltage-drop between two phases with active ion concentrations x₁ and x₂ is Equation (iv):

(iv)

where q is the charge of the active ion, F is the Faraday constant and ΔG_rxn is the change in Gibbs free energy between phases. In practice, the change in Gibbs free energy in Equation (iv) is approximated by the change in the DFT total energy, under the assumption that entropic contributions will have a minimal effect on the free energy differences between phases during cycling.

When studying a phase diagram computationally there are a finite number of phases on the tie-line, thus the profile will not be a continuous smooth line, but a sequence of two-phase regions with constant average voltages. Although the profile will not have the same characteristic curve as an experimental voltage profile, it is still possible to calculate quantities of interest such as theoretical capacity, which is calculated from the maximal difference in active ion concentration between the predicted stable phases. Similarly, the energy density of an electrode is found by integrating the voltage profile between the two endpoint phases.

Beyond calculating the voltage profile, one may further validate a crystal structure against experiment by using DFT to predict its spectroscopic signatures. Many spectroscopic methods, including XAS, EELS (31) and Raman spectroscopy (32), can be readily calculated using DFT to aid characterisation.

Solid-state nuclear magnetic resonance (ssNMR) spectroscopy is a tool for investigating the element-specific local structure of materials, even for the disordered and dynamic systems present in battery materials (33). Due to the complex structures and processes that arise during battery cycling, the usefulness of NMR spectroscopy can be greatly enhanced by applying complementary techniques to aid the assignment of spectra to the local environment of each nucleus. Theoretical methods in DFT are sufficiently mature that the calculation of chemical shielding tensors across a diverse range of inorganic systems is now routine (34).

NMR spectroscopy involves the precise measurement of the response of nuclei in an applied magnetic field to weak oscillating perturbations; for a given pulse scheme, the frequency of perturbing oscillations is adjusted until resonance is achieved, at which point a signal is observed. The frequency of this resonance is a cumulative measure of several competing interactions between the spin of the nucleus and its local environment and, when referenced against a model nucleus, is referred to as the chemical shift. The observed chemical shift in most materials is determined by the nuclear spin interacting with the orbital angular momentum of paired electrons. In Figure 3, such a shift is given for the phases of Li-P which form during cycling of a Li-ion battery with a phosphorus anode (22). The ³¹P chemical shift of each Li_x P_y phase is distinct, as shown by the coloured peaks in the figure for each compound.

Fig. 3.

Whilst the theory for computing magnetic shielding for isolated systems (such as molecules and clusters) was developed in the 1960s and 1970s in the context of quantum chemistry (35), these methods were not easily extendable to solids (36). For periodic systems, such as battery anodes and cathodes, most modern implementations of theoretical ssNMR use DFT and the gauge including projector augmented wave (GIPAW) approach (37–39). It is not only possible to compute the full chemical shielding tensor, but also several other effects that can modify the lineshape of the NMR signal, namely quadrupolar coupling (for spin I >1/2 nuclei), dipolar coupling (which can be simulated directly from the geometry using for example the SIMPSON software package (40)) and J-coupling (interaction of electron spins which can probe chemical bonds directly) (41).

Finally, beyond just characterising the static crystal structure of a battery material, it is also possible to predict the dynamics of ions moving through the material, which is especially useful when studying ionic transport in electrodes and solid electrolytes. The charge and discharge rates are key performance factors in battery design, defining the time required to fully charge a battery and the amount of power it can deliver, respectively. Rate capability is determined by the speed with which the charge carriers can move through the materials. Since both ions and electrons move in a battery, the rate capability depends on both the electronic and ionic conductivity of the materials. While the electrodes in batteries must be mixed electronic-ionic conductors, the electrolyte must be electronically insulating. First principles methods, such as DFT, can be used to study both electronic conductivity and ionic conductivity of battery materials. Electronic conductivity can be assessed from electronic structure calculations (42–44), while ionic conductivity can be calculated using ab initio molecular dynamics (AIMD) or the nudged elastic band (NEB) method (45), as outlined below.

The bulk ionic conductivity, σ(T), of a solid electrolyte can be related to diffusion coefficients via the Nernst-Einstein relation (46) defined as Equation (v):

(v)

where n is the diffusing particle density, e the elementary electron charge, z the ionic charge, k_B the Boltzmann constant, T the temperature, D(T) the ionic diffusivity and H_R the Haven ratio accounting for the correlated ionic motion.

3.3.1 Ab Initio Molecular Dynamics Simulations

One way to compute ionic diffusivity of a given material, using AIMD simulations, combines the first principles aspects of DFT with the ability of molecular dynamics (MD) to model ionic forces and trajectories. Methods to screen the mobility of ions along an MD trajectory include mean square displacement (MSD), mean jump rate (MJR) (47, 48), velocity autocorrelation function (VACF) (49–52) van Hove correlation function (53, 54) and others (55, 56). MSD is the most straightforward and robust and thus the commonly used definition of diffusivity.

One can extract the diffusion coefficient D(T) from the gradient of the MSD, given a well-converged MD trajectory such that the MSD is a linear function of time. Here, the slope of the line of best fit gives the diffusion coefficient D, times twice the dimensionality d of the diffusion (2d * D). For ionic diffusion in three dimensions, d = 3. Depending on the level of mobility of ions in the system, good convergence of the MSD of ions may require long trajectories, for example 50–100 ps, thereby requiring tens of thousands of time steps. As each step involves DFT energy or force evaluations, AIMD can be a computationally demanding process. Two common solutions to this are: (a) to analyse trajectories obtained at elevated temperatures (500–2000 K) to foster higher mobility and faster convergence of the MSD; or (b) to utilise parameterised atomic force-fields to allow faster evaluation of the interatomic forces in the system compared to ab initio methods like DFT. A drawback of parameterised force-fields is non-transferability, so one needs a new set of fitted parameters for the specific set of atoms in a new system.

The activation energy (E_a) for the ionic transport in a given electrolyte or electrode can be obtained from AIMD simulations using the Arrhenius law, Equation (vi):

(vi)

where D₀ is the theoretical maximum diffusivity at infinite temperature, under the assumption that the diffusion mechanism is not temperature dependent and no phase transition occurs. Analysis of the trajectories from the AIMD simulations can also provide useful information on the crystallographic sites with higher occupation probability, while also revealing the preferred ionic conduction pathways between these sites (47, 57, 58).

3.3.2 Nudged Elastic Band Method

Another way to obtain ionic diffusivity from first principles is with optimisation-based methods, through the exploration of minimum energy paths (MEP) describing a set of predefined ionic migration pathways. To this end, the NEB algorithm is often used. Other approaches are also available for transition-state searches, for example the dimer (59), Lanczos (60) and eigenvector-following (EF) (61) methods as well as others (62, 63). Specifically, the NEB method computes the MEP (at 0 K) for a predefined route connecting the initial and final states of the motion of a single ion or a few, concertedly diffusing ions (45, 64). The ion-transport path is divided into intermediate steps (called NEB images), defined by the interpolation of these two end-point states. The NEB images are concurrently optimised by introducing a set of imaginary spring-forces to ensure the harmonic coupling of the consecutive images and a continuous path on the corresponding high-dimensional potential energy surface. Using the climbing-image NEB that maximises the energy of the saddle point(s) on the MEP, one can also locate the transition states, from which activation energies (E_a) are calculated.

In solids, the change of entropy during ionic diffusion is usually negligible and thus activation free energies are typically approximated by their 0 K values. The diffusion rate can then be related to the ionic diffusivity in the dilute carrier limit (65) (i.e. diffusion carriers do not interact) using Equation (vii):

(vii)

where λ is the hop distance between two adjacent sites, g is a geometric factor that depends on the symmetry of the sublattice of interstitial sites, f is the correlation factor, X_D is the concentration of the diffusion-mediating defects, v* is the entropy difference between the initial and final states, the activation energy E_a is the energy difference between the initial and final states, k_B is Boltzmann’s constant and T is the temperature of the simulation.

Static methods, such as NEB, provide computational efficiency over AIMD: NEB requires only a few hundred DFT steps to converge and is accurate within the regime in which the electronic structure of the model system does not change with the ionic migration (66). NEB calculations also allow for quantitative comparison of different migration routes. Nevertheless, NEB is less likely to reveal new conduction mechanisms compared to AIMD, and the complex cooperative conduction mechanisms may not be as straightforward to sample with NEB as with AIMD. Moreover, NEB usually operates in the dilute regime (Equation (vii)), where vacancy defects are manually introduced in the sublattice of the diffusing ions to have a low diffusion carrier concentration and mediate the ionic motion. These artificial defects not only decrease the accuracy of the simulation models, but also impede the integration of the NEB method in high throughput approaches. AIMD, in contrast, would in principle work with any concentration of diffusing ions by readily addressing the self-diffusion limit (67, 68). Given these tradeoffs, a common practice in the literature is therefore to combine AIMD with NEB calculations, specifically by identifying the potential conduction pathways from relatively shorter AIMD trajectories at a selected, elevated temperature and to probe the MEPs to get E_a and compute the other properties relevant to the ionic transport (57, 58, 69–71).

In many cases, as in the high-voltage high-capacity anode material TiNb₂O₇ (TNO), both ionic and electronic conductivities are relevant to the performance of the battery material (72). In this case, density of states (DOS) calculations were used to determine that the electron-doped TNO is metallic, as compared to the pristine TNO. Additional localised electronic states were confirmed in AIMD as a result of bond distortions, thus exemplifying the need in this case for both AIMD and DOS calculations.

Graphite is ubiquitous in contemporary commercial Li-ion batteries. However, alternative anode materials are a highly researched topic, due to graphite’s low capacity (372 mAh g^–1) and tendency for Li plating and subsequent dangerous short-circuiting due to its low operating voltage (76). These factors make graphite anodes unattractive for applications that require high performance and capacity, such as electric vehicles.

Here we highlight developments in predicting high capacity conversion and alloying anodes to replace graphite, based on tin (990 mAh g^–1) and antimony (660 mAh g^–1). Such conversion and alloying anodes undergo a succession of reversible phase transformations during charging and results in their observed higher capacity retention than other conversion and alloying anodes (77).

Both Sn and Sb were previously employed as anodes in Li-ion batteries, showing evidence of conversion reactions, with unknown phases of Li_x Sn and Li_x Sb forming during Li insertion. An AIRSS search for the thermodynamically stable phases of both Li_x Sn and Li_x Sb was conducted (23, 78) in order to understand the voltage profiles and reaction mechanisms of these two alloying anodes. In this case, a new phase Li₂Sn was identified by AIRSS to lie near the convex hull. The resulting voltage profile is compared with experimental measurements in Figure 4.

Fig. 4.

During the cycling process in conversion anodes such as Sn, the material at the anode undergoes several conversion reactions as Li is inserted (77). In the voltage profile shown in Figure 4, the black line is constructed from the ground state phases in the Li-Sn system, which were predicted using AIRSS (23, 78). Each plateau in Figure 4 represents a two-phase region between one ground state Li-Sn alloy and another, until a critical point is reached at which there is a phase transformation (a vertical line) to the next Li-Sn alloy.

The DFT predictions lie within the voltage range of the experiment and are an accurate match to both sets of experimental data by Wang et al. (79). In many cases, the experimental data has less-sharp distinctions between separated phases, due to reactions which appear to occur gradually rather than at a well-defined stoichiometry.

The Li-Sb phase diagram was found to be somewhat simpler, with only two stable phases predicted during cycling: Li₂Sb and Li₃Sb. Two competing polymorphs of Li₃Sb were found and NMR calculations were performed on both to provide a signature of each phase to aid the interpretation of future experiments.

This work on Li-Sn and Li-Sb anodes provided theoretical confirmation of experimental binary phases in this family of conversion anodes and allowed for more concrete evidence of the specific mechanism of Li insertion into these anodes. Furthermore, this study confirmed the new phase of Li₂Sn.

The electrolyte in a battery forms a conductive bridge between the anode and cathode which allows ions to move from one electrode to the other without permitting the flow of electrons. Conventional Li-ion battery architectures use a liquid electrolyte consisting of a Li salt mixture dissolved in an organic solvent. Two prominent safety concerns arise from the use of organic liquid electrolytes (80, 81). The first is that the organic solvent component tends to be flammable and poses a fire hazard when exposed to air if the battery casing is breached (82). The second is that Li dendrites (83, 84) form, which can eventually bridge the gap between the anode and cathode resulting in short-circuiting.

All solid-state batteries attempt to solve these safety issues by replacing the organic electrolyte solutions with solid equivalents, which exhibit high mechanical strength, suppressing dendrite formation, thus enabling the use of the high energy density Li-metal anodes (85, 86). Most proposed solid electrolytes have sufficient mechanical strength, as demonstrated by high throughput screening based on machine learning methods (87).

A key challenge in developing solid electrolytes is finding solids with room temperature (RT) ionic conductivities that approach those of their liquid counterparts. Among several solid electrolyte families identified to date, the thiophosphide ceramics, for example Li₂S-P₂S₅, chemically-doped sulfides, like Li₁₀GeP₂S₁₂ (LGPS) (88) and Li_9.54Si_1.74P_1.44S_11.7Cl_0.3 (89), are known to deliver the highest RT Li-ion conductivities (1.2–2.5 × 10^–2 S cm^–1). Sulfides, however, have high moisture sensitivity and their chemical stability against common electrodes is low, thus limiting their practical use (90). By contrast, oxides like garnets (for example, Li_x La₃M₂O₁₂, where M = zirconium, niobium, tantalum) display notably higher chemical stability than sulfides but exhibit lower ionic conductivities (91). The latter limitation can be partly remedied by a chemical doping with diverse metals, including aluminium, gallium and scandium (92).

High throughput CSP is useful for exploring new superior electrolytes with combined high conductivity and chemical stability. Various studies have performed extensive screening of superionic conductors within databases such as the Materials Project (8), searching for phases with good phase stability, high Li⁺ conductivity, wide band gap and good electrochemical stability (12, 53, 93–95).

Various LGPS-derived compositions were predicted using ab initio calculations through elemental swapping (95), such as Li₁₀(Sn/Si)PS₁₂ and then verified by experimental synthesis and measurements (96, 97). LiAlSO was discovered solely through structure prediction and proposed to be a superionic conductor with AlS₂O₂ layers, which facilitate faster movement of Li-ions, low activation barriers and a wider electrochemical window (94). Similarly, Fujimura et al. (98) presented a high throughput (HT) screening of the chemical phase space for Li_3.5Zn_0.25GeO₄ (LISICON)-type electrolytes. The authors proposed new electrolytes with higher conductivities than the parent LISICON material. Later, Zhu et al. (93) reported a HT screening of the Li-P-S ternary and Li-M-P-S (where M is a non-redox-active element) quaternary chemical spaces and identified two Li superionic conductors, Li₃Y(PS₄)₂ and Li₅PS₄Cl₂. Particularly, Li₃Y(PS₄)₂ is predicted to exhibit a room-temperature Li⁺ conductivity of 2.16 mS cm^–1, which can be further enhanced with aliovalent doping (93). However, these materials are yet to be synthesised.

Following the structure prediction of these new solid electrolyte phases, it is then desirable to use NEB and AIMD simulations to investigate the atomistic origins of their ionic conductivity. For instance, Li-ion transport was elucidated in the sulfide-based electrolytes, Li₇P₃S₁₁ (99), argyrodite Li₆PS₅Cl (48, 53), LGPS (57, 100), Li-Sn-S/Li-Sn-Se (101, 102) and Li-As-S/Li-As-Se alloys (103), Li₃PS₄ (48, 104, 105), Li₄GeS₄ (57, 103) as well as oxides, for example LLZO (71, 106–108), LiTaSiO₅, LiAlSiO₄ (71), Li₄SiO₄−Li₃PO₄ solid mixtures (109) and several others. The problem of identifying solid electrolyte candidates for all solid-state batteries which are air stable and highly conducting can be solved using a combination of structure prediction techniques and atomistic modelling such as AIMD and NEB.

So far, the battery materials we have discussed (Sections 4.1 and 4.2) are based on Li-ion chemistry. However, cost and sustainability are driving research efforts into ‘beyond Li-ion’ batteries. The philosophy presented in Section 3, using CSP and DFT, is straightforward to extend to ‘beyond Li-ion’ chemistries. A prominent example is Na-ion batteries, where Li is replaced with the more earth-abundant Na.

Unlike in Li-ion batteries, graphite shows poor capacity for Na, although other carbonaceous materials offer some promise (110). As such, the success of future Na-ion batteries will rely on the discovery of new anode materials. There are many classes of anode materials which are applicable to Na-ion batteries including two-dimensional transition metal carbides (111) and group V elements (P, As, Sb) (112). Although this review focuses on one Na-ion anode material in particular, structure prediction has been used to predict phases of each of the anode materials in several cases (113–115). In particular, black P shows a high theoretical capacity for Na of 2596 mAh g^–1, corresponding to the formation of Na₃P (22). Here, P acts as an alloying electrode, so its cycling is expected to involve multiple phase transformations. For these reasons, there has been recent focus on understanding sodiation processes in P.

Applying a combination of AIRSS, data mining (22) and a genetic algorithm (25), the convex hull of the Na-P system has been mapped out and is shown in Figure 5. The Na-P system contains a number of stable crystalline phases (coloured black circles in Figure 5) with compositions varying from NaP₇ through Na₃P, and the voltage curve derived from these phases shows good agreement with experimental measurements (25). In addition to these stable phases, there are metastable phases lying close to the convex hull across a range of compositions.

Fig. 5.

By following the structures which fall on or near the convex hull in Figure 5, from least sodiated (pure P) to most sodiated (Na₃P), the calculations predicted many changes in local structure: the layered black P is broken upon successive Na insertion, forming P chains and helices, then dumbbells, which eventually break apart to form isolated P atoms. These structural motifs are distinctive and have characteristic NMR signatures, which can be accurately modelled. In order to confirm this explicitly, ex situ ³¹P solid-state NMR measurements were taken at different points during both the sodiation and desodiation cycle (25). Since contemporary NMR calculations lack a rigorous treatment of paramagnetic contributions to the isotropic shifts, the chemical shift anisotropies were computed for the thermodynamically accessible range of predicted structures to provide a set of chemical environments to screen against experimental measurements. During the reverse cycle when Na is removed from the system, P helices re-formed in a tangled fashion and the original crystalline P was not recovered. Amorphous phases were encountered experimentally on desodiation and, while modelling of amorphous materials is challenging, the local structural features of predicted metastable phases were discovered to be present even in the amorphous structures.

Aside from P, Sn also shows promise as a Na- ion battery anode. Sn presents a lower theoretical capacity for Na (847 mAh g^–1) but offers better capacity retention than P (24). The results of an AIRSS search for Na-Sn phases (24), predicted that insertion of Na into Sn would result in hexagonally layered structures NaSn₃ and NaSn₂, before passing through an amorphous phase of approximate composition Na_1.2Sn, after which a solid-solution consisting of Sn dumbbells surrounded by Na ions would form. The final product, Na₁₅Sn₄, contains isolated Sn atoms surrounded by Na. Importantly, the computational workflow used to study Li and Na-ion batteries is the same and is equally as applicable to conversion anodes for other chemistries.