Computational understanding of Li-ion batteries

Over the last two decades, computational methods have made tremendous advances, and today many key properties of lithium-ion batteries can be accurately predicted by first principles calculations. For this reason, computations have become a cornerstone of battery-related research by providing insight into fundamental processes that are not otherwise accessible, such as ionic diffusion mechanisms and electronic structure effects, as well as a quantitative comparison with experimental results. The aim of this review is to provide an overview of state-of-the-art ab initio approaches for the modelling of battery materials. We consider techniques for the computation of equilibrium cell voltages, 0-Kelvin and finite-temperature voltage profiles, ionic mobility and thermal and electrolyte stability. The strengths and weaknesses of different electronic structure methods, such as DFT+U and hybrid functionals, are discussed in the context of voltage and phase diagram predictions, and we review the merits of lattice models for the evaluation of finite-temperature thermodynamics and kinetics. With such a complete set of methods at hand, first principles calculations of ordered, crystalline solids, i.e., of most electrode materials and solid electrolytes, have become reliable and quantitative. However, the description of molecular materials and disordered or amorphous phases remains an important challenge. We highlight recent exciting progress in this area, especially regarding the modelling of organic electrolytes and solid–electrolyte interfaces.


INTRODUCTION
During the last two decades, lithium-ion battery technology has made possible impressive advances in mobile consumer electronics and electric vehicles. [1][2][3][4] Electrochemical technology for grid-level energy storage additionally bears the potential to enable the transition away from fossil energy towards renewable energy sources. 5,6 To satisfy the increasing demand for low-weight and low-volume batteries with high-energy storage capacity, researchers have been exploring options to increase the specific energy and energy density of lithium-ion batteries. 7 At the same time, the ability to quickly charge and discharge a battery, i.e., the rate capability, is not only critical for the usability of portable devices, such as smart phones, but is also an important factor to render electric vehicles competitive (imagine taking gas would take hours instead of minutes). Despite the need for better performing lithium batteries, the recent battery fires on commercial airplanes are a reminder that safety issues must not be neglected when developing new materials.
Lithium batteries are collections of electrochemical cells, each composed of two electrodes that are separated by an electrolyte. The battery functions by shuttling lithium ions between the electrodes, which typically are intercalation materials, through the electrolyte. The driving force for this process is the difference of the lithium chemical potential in the two electrodes. During discharge, the cathode material (low lithium chemical potential) is electrochemically reduced by intercalating lithium ions from the electrolyte and taking up electrons from an external circuit. Simultaneously, the anode material (high lithium chemical potential) is oxidised. The resulting electric current through the external circuit can be used to perform work, i.e., to run an electronic device. The above process is essentially reversed when the battery is recharged by applying an external potential. To sustain the cell reactions with minimal overpotential, electrode materials have to be good electronic and lithium ionic conductors. The electrolyte, on the other hand, needs to conduct lithium ions but has to be electronically insulating to prevent short circuiting. Today, most commercial batteries employ carbon-based anode materials, though research is in progress to investigate alternatives. [8][9][10] Typical cathode materials are based on rocksalttype lithium transition-metal oxides that provide high-energy densities, such as LiCoO 2 , or polyanionic materials with high-rate capability, such as LiFePO 4 . Commercial electrolytes are organic solutions of lithium salts, 11 but solid lithium ion conductors are under scientific investigation. 12,13 In the past, the development of new materials took place exclusively in the experimental laboratories, often by trial and error, and largely depending on the researchers' experience and intuition. With today's computational capabilities, however, computer simulations have become an integral part of materials design. Atomistic simulations based on first principles, i.e., simulations that are directly based on physical laws, are especially invaluable, as they provide insight into processes on the atomic and electronic scales without requiring any input from experiment, thereby aiding in the interpretation of experimental observations. For applications in materials science, the electronic density-functional theory (DFT) developed by Kohn and coworkers [14][15][16] has been particularly successful.
In recent years, computer-assisted design of entirely new materials has made tremendous progress, not least owing to automated high-throughput calculations that facilitate the rapid screening of thousands of target materials. 17 Computational investigations of specific battery materials have recently been reviewed by Meng and Arroyo-de Dompablo 18 and by Islam and Fisher. 19 Here we will instead focus on the present state of firstprinciples-based methods for the simulation of lithium-ion battery properties. In the following sections, we will review computational approaches to key properties of lithium-ion batteries, namely the calculation of equilibrium voltages and voltage profiles, ionic mobilities and thermal as well as electrochemical stability. Past research efforts in the field have been predominantly geared towards the discovery and optimisation of cathode materials, and thus most examples from the literature cited in the following are for applications to cathode materials. We stress, however, that the techniques are directly transferable to other battery components, such as anode materials, 20-23 solid electrolytes, 12,13,24 and electrode surface coatings. 25

EQUILIBRIUM VOLTAGE
One way to increase the specific energy and energy density of a battery is by increasing the cell voltage as much as possible without exceeding the stability window of the electrolyte. In the following section 'Equilibrium cell voltage from first principles' we will discuss the general first-principles framework for the computation of equilibrium cell voltages. The accuracy of computational voltage predictions and strategies for its improvement are considered in section 'Self-interaction and the accuracy of first-principles voltages'.
Equilibrium cell voltage from first principles The equilibrium lithium intercalation voltage is determined by the difference in lithium chemical potential, μ Li , between cathode and anode where z is the charge that is transferred, and F is the Faraday constant. The lithium chemical potential is the change of the free energy of the electrode material with lithium content. 26,27 Integrating Equation (1) over a finite amount of reaction gives the average voltage as function of the free energy change of the combined anode/cathode reaction (Nernst equation) At low temperatures, the entropic contributions to ΔG r are small, and the reaction free energy can be approximated by the internal energy, ΔG r ≈ΔE r . Within this approximation, the equilibrium voltage of a lithium transition-metal oxide intercalation cathode with composition LiMO 2 and a lithium metal anode with the cell reaction can thus be computed as where the internal energies of the lithiated and delithiated phases, E(Li x1 MO 2 ) and E(Li x2 MO 2 ), and of metallic (body-centered cubic) lithium, E(Li), can be obtained from first principles. Obviously, this approximation is not limited to intercalation electrodes, but can also be applied to conversion or displacement reactions.
At this point it is worth to step back and look at the first principles approximation to the equilibrium voltage, Equation (4). All that is required to compute the voltage are three independent first principles calculations for Li x1 MO 2 , Li x2 MO 2 , and Li, and the energy of BCC lithium is independent of the cathode material and hence only needs to be computed once. This means, the average intercalation voltage of, e.g., layered LiCoO 2 can be estimated simply based on the results of DFT calculations of LiCoO 2 and the delithiated CoO 2 , calculations that can be done within minutes on a current desktop computer. Note that if phases with intermediate lithium concentrations are known, Equation (4) can be used to compute a piece-wise approximation to the voltage curve. However, the real challenge lies in the determination of the relevant thermodynamically stable phases Li x1 MO 2 and Li x2 MO 2 and their respective crystal structures, a topic that will be addressed in section 'Voltage profiles'.
Self-interaction and the accuracy of first-principles voltages Most early DFT calculations for battery materials were based either on the local density approximation (LDA) 15 or the generalised gradient approximation (GGA). 28 Aydinol et al. 27,29,30 (for the case of LDA) and Deiss et al. 31 (GGA) demonstrated that the approach outlined above reproduces experimental voltage trends. However, a systematic underestimation of experimentally measured voltages of layered lithium transition-metal oxides by up to 1.0 V was found. Despite this large systematic error, DFT calculations were instrumental for the understanding of the fundamental electronic structure origin of redox levels. To this extent, early LDA calculations revealed that the average intercalation voltage of layered LiMO 2 (M = first-row transition metals) increases with the atomic number of M, as a result of the increasingly covalent character of the M-O bond in these materials. 27,32 The same general trend was later also confirmed for polyoxianionic cathode materials. 33,34 In addition, LDA predicted a decrease in the intercalation voltage with the electronegativity of the anion, V(LiMO 2 ) 4V(LiMS 2 ) 4V(LiMSe 2 ). 27,35 The impact of the host structure on the voltage is generally found to be small in comparison to the voltage variations upon cation doping and exchange of the transition-metal species. 27,35 The reason for the failure of DFT to predict the voltages of lithium transition-metal oxides quantitatively correct can be traced back to a general problem of (semi-)local DFT. The meanfield approximation to the electrostatic interaction of the electrons in DFT introduces a spurious self-interaction, i.e., the interaction of each electron with itself that is not fully canceled out by LDA and GGA functionals. This self-interaction error results in an artificial delocalisation of the electrons that leads to significant errors in systems with strongly correlated electrons. 36 The true electron density of many transition-metal oxides, for example, exhibits strong features of the individual transition-metal atoms with localised d electrons at the metal centers that are not captured by conventional DFT calculations. In fact, GGA incorrectly predicts metallic ground states for many insulating transition-metal oxides that are relevant for lithium ion batteries. 36 Applications that solely rely on relative energies, for example the energetic comparison of different polymorphs, may benefit from error cancellation which effectively reduces the self-interaction error. However, modelling redox reactions requires computing the energy difference of materials with very different electronic configurations, e.g., transition-metal oxides in their oxidised and reduced states, and thus no such error cancellation can be hoped for.
A practical approach to correct the self-interaction error is offered by the DFT+U method, 37-39 which essentially introduces a penalty term for partial occupations, favouring the disproportionation into fully occupied and empty states. 36 The method lends ideas from the Hubbard model Hamiltonian that was developed to investigate the electronic structure of strongly correlated materials. The magnitude of the DFT+U correction is controlled by a system and implementation specific Hubbard U parameter that can either be determined self-consistently from linear response theory 36,40 or by fit to reference band structures or formation energies. 41,42 For a given Hubbard U, DFT+U calculations are no more computationally demanding than conventional DFT calculations.
Zhou et al. 41 showed that DFT+U calculations using selfconsistent U parameters, significantly improve the accuracy of predicted voltages over pure GGA for transition-metal oxides and polyanionic materials, reducing the deviation from the experimental voltages to about ± 0.1 V. The authors report that self-consistent U parameters determined for either the oxidised or the reduced cathode materials result in comparable accuracy. In addition, DFT+U also correctly describes charge localisation and disproportionation, which is not captured by pure GGA/LDA calculations. 43 Yahia et al. assessed the impact of different U values on the relative phase stability of different LiMSO 4 F (M = Fe, Mn) polymorphs and found only a small dependence, demonstrating the transferability of the U parameter. 44 Another avenue to reduce the self-interaction error of LDA/GGA is to exactly evaluate the exchange energy, i.e., the non-classical contribution to the electron-electron interaction, based on the expressions from the Hartree-Fock (HF) method. 45,46 Exchangecorrelation functionals that implement this approach are called hybrid functionals, as they essentially are a blend of DFT and HF. Instead of requiring additional parameters for selected atomic orbitals (as in DFT+U calculations), the only adjustable parameter in hybrid functionals is the amount of exact exchange, which is only weakly system dependent. 46,47 Note that the long-ranged nature of the electrostatic potential makes it computationally demanding to evaluate the exact exchange energy in periodic calculations. Range-limited hybrid functionals, such as the Heyd-Scuseria-Ernzerhof (HSE) functional, 48,49 reduce the computational effort to some extent by considering only short-ranged contributions. Nevertheless, hybrid-functional calculations scale less favourable with the system size as conventional DFT calculations and are therefore computationally significantly more demanding.
Chevrier et al. 50 showed that hybrid functional DFT calculations using the HSE functional predict redox potentials with a similar accuracy as GGA+U, albeit without the need for additional parameters. Even with HSE significant discrepancies can exist between predicted and measured voltage for systems with complex electronic structure, such as Li x CoO 2 , which undergoes a metal-insulator transition when delithiated. 51,52 Seo et al. 47 found that system-specific adjustment of the amount of exact exchange against reference band gaps from spectroscopy can further improve the accuracy of equilibrium voltages and voltage profiles computed with hybrid functional DFT. A comparison of equilibrium voltages of lithium transition-metal oxides and phosphates computed with GGA, GGA+U and HSE06 are shown in Figure 1.

VOLTAGE PROFILES
Many lithium-cathode materials exhibit stable phases at intermediate lithium concentrations, such as lithium-vacancy orderings in intercalation materials or different atomic orderings in alloys. If the structures of all phases are known, then the equilibrium voltage between the intermediate phases can be computed using Equation (4) of the previous section.
An early example is the calculation of the 0-K voltage profile of the Li-Sn alloy by Courtney et al. 53

0-K intercalation voltage curves
The relevant quantity to compare the stability of different phases at 0 K is the formation energy with respect to stable reference materials. For the example of an intercalation voltage curve for a lithium transition-metal oxide with formula unit LiMO 2 , the formation energy of any structure with intermediate lithium content can be expressed as where E is the internal (DFT) energy and the fully lithiated LiMO 2 and delithiated MO 2 phases are the relative energy reference.
The formation energies of all Li x MO 2 phases that are (at 0 K) thermodynamically stable compared with the reference phases lie on the lower convex hull of E f versus composition x. 55,56 Once the convex hull construction is available, a piecewise voltage profile can be obtained by evaluating the equilibrium voltage, Equation (4), between phases of adjacent lithium concentrations. Figure 2a shows an example of a formation energy hull construction and the corresponding 0 K voltage curve for sodium intercalation in Na x MnO 2 , a system with a particularly large number of stable intermediate phases. 57 Although technically not a lithium battery system, Figure 2a perfectly illustrates the relationship between the formation-energy hull and the 0-K voltage profile, and it includes the measured voltage profile for reference. As seen in the figure, the computational voltage profile (red solid lines) traces the experimental profile (black dotted lines) well, indicating that phases that are stable at 0 K also dominate at operation conditions. Given a set of trial atomic structures, the convex hull construction can be used to identify the stable phases among those structures, but the trial structures need to be obtained from somewhere in the first place. In the case of conversion materials, likely candidates for intermediate structures can sometimes be guessed based on the structures of known phases and chemical intuition. 58 However, for intercalation systems, all intermediate phases belong to the same host structure, i.e., lithium ions and vacancies occupy a common sublattice. Although this general problem of finding the distribution of a species over available lattice sites can be solved by the cluster expansion technique, as has been demonstrated for Li x CoO 2 (refs 55,59) and Li x NiO 2 , 60 Figure 1. Comparison of the errors in intercalation voltages calculated with GGA (red circles), GGA+U (green triangles), and the HSE hybrid functional (blue squares). 50 The Hubbard U value for Ti is equal to zero, so that GGA = GGA+U. (Copyright: American Physical Society).
one usually defers to a simpler enumeration of likely arrangements using a technique that was originally developed by Hart et al. [61][62][63] for the enumeration of atomic orderings in alloys. However, the low-vacancy and low-lithium-concentration regimes require large simulations cells, and for intermediate concentrations the number of possible lithium/vacancy orderings even in small cells may become exceedingly large due to the combinatorial explosion, preventing the first principles energy evaluation of all configurations. Intuitively, one could argue that the most stable phases should be those lithium/vacancy orderings that distribute both species homogeneously over the available sites, as such arrangements minimise the electrostatic repulsion between positively charged lithium ions. Following this line of reasoning, Hautier et al. proposed to rank the enumerated configurations by their electrostatic energy, assuming fixed atomic charges, before considering only the N lowest energy configurations for first principles calculations. 64 This approach is generally useful to determine atomic orderings for crystal structures with fractional site occupancies: For example, Kim et al. employed the enumeration technique in conjunction with filtering by electrostatic energy to identify stable LiMn 0.5 Fe 0.375 Mg 0.125 BO 3 phases in which Mn, Fe and Mg share a common sublattice. 65 We note in passing that the Python Materials Genomics package provides a convenient interface for the enumeration of atomic configurations and allows filtering by electrostatic energy. 66 Finite-temperature voltage profiles Much of the discrepancies between the experimental voltage curve and the piecewise computational approximation in Figure 2a arises from temperature effects that were neglected in the previous section. Indeed, some of the sodium/vacancy orderings (shown as black crosses in Figure 2a) that are not ground states are within a small energy above the convex hull, so that they may become accessible at finite temperatures.
Including temperature effects in computational voltage profiles requires knowledge of the internal energy and the entropy of the thermalised system to minimise the temperature-dependent free energy of the system. A standard technique for the sampling of the configurational space at some given conditions is the Metropolis Monte-Carlo (MC) method, which stochastically samples the system's states with their correct (Boltzmann) probabilities in a set statistical ensemble. 67,68 However, for the convergence of thermodynamic ensemble averages, MC simulations typically require on the order of millions of energy evaluations and atomic configurations that are too large to be calculated easily by DFT.
The cluster expansion [69][70][71] previously mentioned is an elegant approach to map accurate DFT energy models onto a simpler Hamiltonian, which captures the dependence of the energy on the Li/vacancy distribution, and can be evaluated fast enough to be used in MC simulations. In the case of lithium/vacancy orderings that only differ in the arrangement of lithium ions and vacancies on their respective sublattice, the largest contribution to entropy can be expected to be due to the configurational degrees of freedom. Although the magnitude of the vibrational entropy may be significant, its variations across phases of the same chemical species are typically small so that it does not much affect the relative phase stability. 72 In addition, electronic contributions to the configurational entropy due to the ordering of transitionmetal ions of the same chemical species but in different oxidation states (electron-hole ordering) may become significant in the presence of highly localised electrons. 73,74 Fortunately, the atomic interactions in crystalline solids can be well discretised to lattice models, which allow rapid energy evaluations of configurations of several thousand atoms. In the cluster expansion method, sometimes also referred to as generalised Ising model or lattice gas model, [69][70][71] the total energy of an atomic configuration s = {σ i }, i.e., a particular lithium/vacancy ordering, is expressed as an expansion over energy contributions by clusters {α} of lattice sites where J 0 is a constant shift and the expansion coefficients J α are the effective cluster interactions (ECIs). For most relevant systems, the ECIs decay rapidly with the number of sites within the cluster (n-body interaction) and with their distance from each other (interaction range), so that the expansion (Equation (6)) quickly converges and may be truncated in good approximation. Several standard methods are available to determine the non-zero ECIs {J α } via fit to a small number of first-principles reference energies. [75][76][77] The cluster expansion method and related lattice-based models have been applied to various battery materials. Early examples from the literature are the simulation of the room-temperature voltage profile of the Li-Al alloy by Reimers and Dahn, 78 and the investigation of lithium-vacancy orderings in Li x CoO 2 by Van der Ven et al. 55,59 Wolverton and Zunger also applied first-principlesbased cluster expansion to cation (Li/Co) and lithium-vacancy ordering in Li x CoO 2 , reproducing the experimentally observed layered ground state and reporting the temperature-dependence of the voltage profile. 79,80 Van der Ven and Ceder used MC simulations based on a cluster expansion Hamiltonian to compute the temperature dependence of the lithium intercalation voltage profile of Li(Ni 1/2 Mn 1/2 )O 2 (Figure 2b). 81 As seen in Figure 2b, temperature effects mainly smooth out the voltage steps that are present in the 0-K approximation, resulting in a continuous voltage slope instead, so that the difference between finitetemperature and 0-K voltage profiles are often small. As such, the 0-K approximation is usually sufficient to gain qualitative insight into the charge mechanism and to estimate the energy density. More recently, Lee and Persson 82 employed a cluster expansion model to determine the structural phases that give rise to the different features in the voltage profile of Li x Ni 0.5 Mn 1.5 O 4 . Yu et al. 83 calculated the phase diagram of spinel Li x Cu y TiS 2 , a material that undergoes lithium displacement reactions, using MC simulations based on a cluster expansion Hamiltonian as part of a multiscale simulation effort.

IONIC MOBILITY
Apart from the energy density, the rate capability is another key factor for the design of new battery materials. The time that is required for a full charge is, for example in the case of a cell phone or electric vehicle, critical for the usability of the device. The discharge rate, on the other hand, determines the amount of power that can be delivered by the battery.
One rate-critical process in lithium intercalation batteries is the extraction of lithium atoms from and their reinsertion into the host structures of the electrode materials. The intercalation rate can either be limited by electric conductivity or ionic conductivity. In this section, we shall focus on the latter, i.e., on the mobility of lithium ions within the intercalation host. Note that the same considerations also apply for ionic diffusion through solid electrolytes.
On the microscopic scale, the chemical diffusivity can be expressed in terms of the Einstein relation 84,85 where d is the dimensionality of the diffusion, e.g., d = 2 for lithium diffusion in layered oxides, R(t) is the total displacement of all N diffusing particles during the time period t, r ! i ðtÞ are the individual particle displacements, and ::: h i indicates the statistical average. The thermodynamic factor Θ accounts for the fluctuations of the chemical potential of the diffusing species, μ, with the concentration and can be expressed as Note that the thermodynamic factor is related to the negative slope of the voltage profile, ∂μ/∂c = − xF ∂V/∂c, as evident by comparison of Equation (8) with the expression of the voltage, Equation (1), in terms of the lithium chemical potential. In the following section 'Direct simulation of the diffusion dynamics', we consider the evaluation of the diffusivity by direct ab initio molecular dynamics simulation of the ionic diffusion. However, in many cases transition-state theory based on 0-K diffusion barriers provides sufficient insight into lithium migration (section 'Diffusivity based on 0-K migration barriers'). In section 'Conceptual insight into lithium diffusion', we discuss how the understanding of lithium diffusion on the atomic scale can be used to devise design criteria to optimise macroscopic lithium transport.
Direct simulation of the diffusion dynamics The migration of a lithium ion from one site to another is an activated process with a free energy barrier. One way to gain insight into lithium ion diffusion and its microscopic diffusion mechanisms is by direct ab initio molecular dynamics (AIMD) simulations. In AIMD simulations, the atomic forces from first principles (i.e., quantum mechanics) methods are used to propagate the atoms in the system according to the laws of classical mechanics. For a general introduction to the subject we refer the reader to standard text books. 86,87 AIMD simulations can most readily address the limit of self diffusion, i.e., the case where no concentration gradient is present, for which the thermodynamic factor, Θ of Equation (8), is equal to one, as the explicit evaluation of Θ in AIMD simulations is not straightforward. In the dilute limit, the calculation of the diffusivity further simplifies, as the time-correlations between individual particle positions can be neglected, r ! i r ! j % 0. The ensemble average of the total displacement in Equation (7) can, hence, be replaced by the atomic mean squared displacement (MSD), r 2 , so that the diffusion coefficient for Θ≈1 becomes This approximation is beneficial, as averaging over all N atoms improves the sampling statistics, and hence the ensemble average converges more rapidly with the simulation time. Note that D * is also called the tracer diffusion coefficient. As pointed out by Alder et al., 88 the equivalent long-time limit of the change of the MSD with time converges faster and hence this expression is usually used in practice. Equation (10) provides an efficient approach to obtain the diffusivity from AIMD simulations at constant temperature and volume, i.e., simulations in the canonical (NVT) statistical ensemble. One simply has to compute the MSD and in a plot of the MSD against the simulation time the slope is 2d D.
Note that for systems with highly correlated diffusion, the tracer diffusion coefficient may not be a good approximation, and the diffusivity should be directly evaluated according to Equation (7), i.e., based on the ensemble average of the total displacement, requiring longer MD trajectories to achieve convergence. An indepth discussion of the different approximations to the diffusivity and their validity and relationship can be found in reference. 89 The direct AIMD simulation of the lithium diffusivity at operation temperature is generally challenging because of the low lithium diffusivities. For the case of lithium diffusion through typical electrode materials, D is of the order of 10 − 10 to 10 − 6 cm 2 /s at 400 K, 90 i.e., the mean displacement an atom experiences in one picosecond (at least 500 AIMD steps) is 0.001 to 0.1 Å. AIMD simulations of several nanoseconds would be required to observe ionic migration, yet longer trajectories to converge the value of the diffusivity. However, if the diffusion mechanism is independent of the temperature, then the diffusivity is often found to follow the Arrhenius law where E a is the activation energy of the diffusion. In practice, the diffusivity can therefore be evaluated at elevated temperatures (500-1500 K) at which shorter AIMD trajectories suffice, and values at lower temperatures can be obtained by extrapolation of log D.
For comparison with experimental results, it is useful to relate the lithium ion diffusivity to the ionic conductivity σ via the Nernst-Einstein relation where N Li is the number of lithium ions, V is the volume of the simulation cell, e is the electronic charge, and T is the temperature.
In the context of lithium ion batteries, Yang and Tse 91 reported AIMD simulations of lithium diffusion in LiFePO 4 , identifying a diffusion mechanism that involves the creation of Li − Fe anti-sites. Mo et al. used AIMD simulations to estimate the lithium diffusivity in Li 10 GeP 2 S 12 , a super ionic conductor material and prospective solid electrolyte, and report a dependence of the diffusivity on the lattice direction, originating from differences in the diffusion pathways. 12 An AIMD trajectory from this work is shown in Figure 3a. The same authors employed AIMD simulations to identify the sodium diffusion pathways in P2-Na x CoO 2 and related materials. 92 Hao and Wolverton investigated lithium transport in the amorphous electrode coating materials Al 2 O 3 and AlF 3 using AIMD simulations. 93 Xiao et al. used AIMD simulations to assess lithium diffusion in the transition-metal layer of lithium rich Li 2 MnO 3 , observing a diffusion mechanism perpendicular to the layer. 94 Diffusivity based on 0 K migration barriers Although AIMD simulations in principle provide a parameter-free route to calculate ionic diffusivities, they are computationally quite demanding. Each step of an AIMD simulation basically is a separate DFT calculation, and the convergence of the diffusivity typically takes on the order of tens of thousands of AIMD steps. 12 Therefore, brute-force AIMD simulations should be the last resort, in cases where no further information about the lithium diffusion mechanism in the system is available, or when the lithium diffusion mechanism is too complex to capture with a simple hopping mechanism.
Consider the microscopic mechanism of lithium hopping, which is at the origin of diffusion. Each such individual hopping event requires a Gibbs free energy of activation, ΔG ‡ , that is given by the free energy difference of the initial state (i.e., the lithium ion in its original site) and the energetically highest state that has to be overcome during the diffusion, the transition state. According to transition-state theory, 95,96 the rate k at which the hopping in which ν* is a temperature-dependent effective attempt frequency. When the hopping distance between adjacent sites, a, is known, the chemical diffusivity in the dilute carrier limit can be approximately obtained from the rate as D(T) = a 2 k(T). 97 However, for systems that follow the Arrhenius law, Equation (10), the pre-exponential coefficient ν * and the activation energy can be assumed to be independent of the temperature. By neglecting the change of entropy during the diffusion, ΔG ‡ can be approximated by a 0-K activation energy ΔE ‡ = E ‡ − E i , where E ‡ and E i are the energies of the transition state and the initial state, respectively. Typical values for the prefactor ν* in Equation (13) are 10 11 to 10 13 s − 1 . 98 On the basis of the knowledge of the microscopic diffusion rates k of Equation (13), the actual ionic diffusivity, Equation (7), can be approximated as D≈g·a 2 ·k, where g is a geometric factor and a is the hopping distance between two adjacent sites. 99 The geometric factor is usually close to 1 and therefore often neglected, 56 however, an explicit calculation of the diffusivity is feasible with kinetic Monte-Carlo (kMC) simulations on lattice models, 85,100,101 such as the ones discussed in section 'Finite temperature voltage profiles'. By combining the energies of transition states with the cluster-expanded energies of lattice configurations, a complete study of lithium diffusion including composition and ordering dependence can be undertaken. 85 Lattice-based Monte-Carlo simulations have the additional advantage over AIMD simulations that the thermodynamic factor, Equation (8), can be directly evaluated in grand-canonical simulations. 98 An efficient algorithm for the computation of transition-state energies is the nudged elastic band (NEB) method. 102,103 The NEB algorithm requires the initial and final states of the diffusion as input, from which it generates a number of intermediate states, the images, by linear interpolation. The minimum energy path (MEP) connecting initial and final states is then determined by concurrent minimisation of the atomic forces in all images subject to a harmonic coupling between neighbouring images. See Figure 3b for example MEPs for lithium diffusion in LiTiS 2 . In our experience, the NEB method is very robust and reliably converges the MEP, as long as the electronic structure of the system does not significantly change during the migration. For practical purposes, this means that GGA+U calculations (section 'Self-interaction and the accuracy of first-principles voltages') tend to be more problematic, as electrons might be localised at different atomic centers along the diffusion path, which may result in the simultaneous diffusion of a polaron that gives rise to an additional charge-transfer barrier. 104 A common strategy to decouple the ionic diffusion barrier from the effects of electron-hole interactions is to turn to plain GGA calculations in which electrons are more delocalised (see also section 'Self-interaction and the accuracy of first-principles voltages'). 105 NEB calculations have clarified the diffusion mechanism in several important cathode materials. As quantitative data and mechanistic information is difficult to obtain experimentally, lithium diffusion is an area where computation has been particularly useful. In early work, Van der Ven and Ceder 90,106 employed NEB calculations to understand the lithium diffusion mechanism in layered LiCoO 2 , identifying a high-barrier diffusion path through an oxygen-oxygen bond and a low-barrier di-vacancy pathway via a tetrahedrally coordinated activated state. This di-vacancy mechanism controls lithium diffusion at all practical concentrations and has become generally accepted as the reason why layered cathodes have good lithium mobility. The authors subsequently carried out kMC simulations at different lithium concentrations, predicting maximal lithium diffusivity in partially delithiated phases in agreement with experiment. Morgan et al. 97 estimated the effect of cation substitution on the lithium diffusivity in olivine LiMPO 4 (M = Mn, Fe, Co, Ni) from NEB diffusion barriers and TST, finding a one-dimensional diffusion mechanism and generally low activation barriers on the order of 100-300 meV, which set the stage for understanding the sizedependence behaviour of lithium diffusion in LiFePO 4 (ref. 107) and led to the evaluation of the very high rate LiFePO 4 . 108 Van der Ven et al. 98 employed MC and kMC simulations to investigate lithium diffusion in Li x TiS 2 accounting for the concentration-dependent thermodynamic factor, Θ. Θ is found to vary four orders of magnitude between the dilute lithium and the dilute vacancy limit due to lithium-lithium interaction. Kang and coworkers employed AIMD simulations at high temperature to identify the lithium diffusion pathways in Al-doped LiGe 2 (PO 4 ) 3 , a solid electrolyte, and subsequently used NEB calculations to converge diffusion barriers and to estimate the lithium diffusivity. 109 The simulations predict that Al doping introduces an alternative diffusion mechanism that enhances the lithium diffusivity compared with the undoped material. Du et al. 110 computed the lithium migration barrier in the crystallographic c direction of Li 10 GeP 2 S 12 with NEB. 12 The result of 230 meV is slightly larger than the value of 170 meV obtained from AIMD by Mo et al. (see previous section). 12,110 Finally, a recent review of lattice-based simulations of lithium diffusion in intercalation materials can be found in ref. 56 Conceptual insight into lithium diffusion The previous two sections dealt with different approaches for the computational estimation of lithium mobility by simulating lithium diffusion. Although these computational tools are invaluable to understand new materials, sometimes useful concepts can be identified for entire materials classes.
Conceptually, lithium diffusion is best understood for rocksalttype lithium transition-metal oxides. In fully lithiated LiMO 2 phases, each cation (lithium or transition metal (TM)) is octahedrally coordinated by six oxygen atoms, and lithium diffusion from one octahedral site to another octahedral site takes place via a tetrahedral activated state (o-t-o migration). 106 Figure 3b depicts the MEPs for lithium o-t-o migration through the layered and spinel structures at different local lithium concentrations, and the local minimum halfway corresponds to the activated state. 56 The lower migration barrier in the layered structure (red curve in the top panel of Figure 3b) corresponds to the di-vacancy mechanism mentioned above, i.e., the additional lithium site adjacent to the activated tetrahedral state is vacant. 90,106 As shown for the specific case of LiTiS 2 in Figure 3b (top), the di-vacancy mechanism reduces the migration barrier by around 300 meV or almost 50%, which makes it the dominant mechanism in this material. Vacancies lower the migration barrier as the energy of the activated state is mainly determined by the electrostatic repulsion between the activated lithium ion and its neighbouring cations. 111 In the layered structure, every diffusion channel passes along a (TM) ion, i.e., each activated tetrahedral lithium atom has one neighbouring site that is occupied by a TM (1-TM diffusion channels). As a consequence of the electrostatic repulsion between the migrating lithium and the TM, the distance to this neighbouring TM ion (as controlled by the lattice parameters) and the valence of the TM species can be tuned to optimise lithium mobility, as demonstrated by Kang et al. 112 for nickel based layered LiMO 2 .
The spinel structure also exhibits diffusion channels without an adjacent TM ion, and as a consequence a tri-vacancy diffusion mechanism without any adjacent cations becomes available upon lithium extraction (0-TM channels). The migration barrier associated with 0-TM channels is yet lower (green line in the bottom panel of Figure 3b) and is mostly independent of the TM Computational understanding of Li-ion batteries A Urban et al species in the material. Lee et al. 113 showed that 0-TM channels are also present in cation-disordered materials, but that an excess of lithium is required so that they form a percolating network throughout the material and can be utilised for fast macroscopic lithium transport. The lithium concentrations that enable percolating 0-TM channels for different degrees of cation mixing in layered LiMO 2 are mapped as the blue area in Figure 3c. As shown in the figure, the percolation threshold is around 10% excess lithium in the disordered structure (100% cation mixing) and has a minimum of about 6% in partially cation-mixed layered materials (50% cation mixing). Urban et al. subsequently generalised this percolation model to other rocksalt-type lithium transition-metal oxides with arbitrary cation order, identifying partially cationdisordered spinel structures as a class of fast lithium ion conductors. 114 THERMAL AND ELECTROCHEMICAL STABILITY Although the motivation for the development of new battery materials is foremost to improve the performance of lithium-ion batteries, the safety of the technology must be maintained. In their charged state, many lithium battery materials tend to be thermally unstable. In particular, oxide-based cathode materials, such as LiMO 2 (M = Mn, Co, Ni), experience, as they become more oxidised, a thermodynamic driving force for electrochemical reduction by release of oxygen gas a reaction that is potentially exothermic and may thus lead to thermal runaway and ignition of the electrolyte. Oxygen release is of particular importance at particle surfaces, where it may result in the formation of other surface phases or at high temperature lead to combustion of the electrolyte. 115 Decomposition will also occur when the cell voltage exceeds the stability window of the electrolyte, resulting in electrochemical oxidation of the electrolyte molecules. Computational handles to estimate the stability of the cathode and the electrolyte under operation conditions are therefore an essential part of in silico design of new battery components.
Thermal stability of cathode materials The kinetics of cathode decomposition and reaction with the electrolyte is complex and outside of the capabilities of today's ab initio methods, though interesting work to understand the reaction paths that lead to anode/electrolyte passivation have been explored with ab initio approaches. 116 Hence, focus has been on understanding better the driving force for thermal decomposition to other phases. The first step in assessing the thermal stability is the computational construction of a phase diagram that captures the more reduced phases well. As discussed in the section '0-K intercalation voltage curves' for the case of voltage profiles, a 0-K phase diagram can then be obtained by constructing the lower convex hull of the formation energies of all relevant phases. Not unlike voltage calculations, the construction of phase diagrams requires accurate free-energy differences of electronically different phases. Hence, it is clear that the spurious DFT self-interaction (section 'Self-interaction and the accuracy of first-principles voltages') will introduce significant errors in phase diagrams of transition-metal oxides due to the artificial delocalisation of the transition-metal d electrons. Another large error in the formation energies of oxides arises from an artificial stabilisation of the O 2 molecule in GGA and LDA calculations. These systematic errors in local DFT need to be addressed as thermal decomposition originates from the competition of the charged state with more reducing phases. Wang et al. analysed these two error contributions and proposed a correction scheme based on (i) a constant energy shift of the DFT O 2 energy to account for overbinding, and (ii) a Hubbard U correction (section 'Self-interaction and the accuracy of first-principles voltages') using U values fitted to experimental oxide formation energies. 42 Following this approach it was possible to reduce the error in the formation energies of transition-metal oxides from up to 1.0 eV with GGA to less than 0.1 eV in most cases. 42 For solid phases, 0 K phase diagrams are generally a good first approximation for the true, finite temperature phase diagrams. For example in the case of LiFePO 4 , Ong et al. found that the 0 K Li-Fe-P-O phase diagram well predicts the experimentally known stable phases. 117 However, temperature dependence has to be reintroduced to estimate thermal stability. In particular, the oxygen release in reaction (14) results in a significant increase in entropy that has to be accounted for, for instance by expressing the Gibbs free energy of the oxygen release reaction (14) as 118 where the entropy for oxygen gas at a given temperature, S p 0 O2 ðTÞ, is obtained from experimental thermochemical data for a reference oxygen partial pressure p 0 , and the energy of the O 2   Figure 4a), 118 finding that delithiated layered phases are generally metastable and that phase separation into the spinel structure LiM 2 O 4 and either the layered LiMO 2 or oxygen-deficient structures is thermodynamically preferred. The reaction enthalpies computed according to Equation (15) were generally within 10-20 meV/formula unit of the experimentally measured value for the case of Li-Ni-O 2 .
The approximation of the reaction free energy in Equation (15) does not acknowledge whether the reaction conditions are oxidising or reducing, as controlled by the gaseous environment and the temperature. To describe phase equilibria with respect to the ambient conditions in a system that is open to oxygen uptake or release, the oxygen grand potential has to be considered 117 where N O2 is the number of oxygen molecules and the oxygen chemical potential, μ O2 , depends on the oxygen partial pressure, p O2 , and the temperature: The oxygen enthalpy, H p 0 O2 , can be approximated using the correction scheme by Wang et al. 42 119 Their computations predict, in agreement with experiment, that the manganese phosphate material reduces at lower temperature than the iron-based material, and the predicted oxygen gas evolution versus temperature is shown in Figure 4b. We note, however, that oxygen evolution at the particle surface may be kinetically limited, as pointed out by Mo et al. 120 for the example of lithium peroxide.
Electrochemical stability of electrolytes Electrolytes in lithium-ion batteries are typically solutions of lithium salts in organic solvents containing additional additives, for example, to enhance the solubility or to create stabilising passivation layers on the electrodes. For the electrolyte to be stable at operation conditions, none of its components may participate in electrochemical reactions with the electrodes. In other words, the cell voltage must at no time leave the voltage range over which the electrolyte is stable, otherwise either reduction or oxidation of the electrolyte molecules would occur.
The redox potentials of the electrolyte molecules are determined by the energy that is required to release and take up electrons, i.e., by the ionisation potential and the electron affinity, both of which can be directly evaluated by first principles calculations of isolated charged molecules. 121 However, the relevant energies are not those of molecules in the gas phase, but instead the solvation energy of all involved species needs to be accounted for. The Gibbs free energies of reduction and oxidation of the solvated molecule M that determine the redox potentials are reduction : where '(s)' indicates species in solution. Zhang et al. 122 proposed an indirect method to compute the electrolyte redox potentials based on thermodynamic cycles. A schematic of an essentially identical approach is shown in Figure 5a. In Zhang's method, the solvation energy of the molecular species is approximated by a continuum solvent model, and the overall reaction free energy in solution is expressed in terms of the free energies in the gas phase and the solvation energies (see equations in Figure 5a). The reference energy of the Li/Li + redox couple can be calculated based on a similar thermodynamic cycle, 122,123 or can be obtained from measurement versus the standard hydrogen electrode. 124,125 Note that the free energy of electrons in the gas phase, G[e − (g)], does not occur in the final expression of the voltage, as it is canceled out when the electrolyte potential is referenced to the Li/Li + redox couple. Shao et al. 123 employed Zhang's approach to the calculation of the electrochemical windows of sulfone-based electrolytes with explicit evaluation of the Li/Li + redox potential in solution, comparing the accuracy of different first-principles methods and continuum solvation models. Regarding the solvation models, the authors report that the polarised continuum model (PCM) 126 results in the smallest errors compared with experiment. The most accurate electronic structure method for the prediction of the electrolyte potentials was found to be Møller-Plesset perturbation theory (MP2) with errors in the range of 0.1-0.5 V, whereas DFT (GGA and hybrid functionals) resulted in errors of up to 1.5 V. 123 The validity of approximating solvation effects using continuum models was further assessed by Ong et al. 127 , who compared the electrochemical stability windows of ionic liquids predicted by a molecule-based estimate using PCM with the results of direct classical molecular dynamics simulations of the liquid phase followed by DFT. The authors find that while the stability is generally overestimated by the PCM approach, the identified trends are nonetheless in agreement with explicit solvent calculations.
Wang et al. 128 applied Zhang's method to calculate the oxidation potentials of redox shuttle additives. The authors find a strong correlation of the oxidation potential with the energy of the highest occupied molecular orbital (HOMO). Cheng et al. used Zhang's method for the high-throughput screening of around 1,400 organic molecules, referencing the free energy differences to the experimental voltage of the Li/Li + redox couple (1.24 V). 124 The comparison of the results for this large number of molecules also confirmed that the energies of the HOMO and lowest unoccupied molecular orbital (LUMO) correlate well with the ionisation potential and the electron affinity (Figure 5b), 121,124 thus eliminating the need to calculate the energy of the reduced and oxidised species in a rapid screening approach.
Borodin et al. investigated the effect of nearby anions on the oxidation stability of solvent molecules and found that hydrogen or fluorine abstraction may significantly reduce the oxidation potentials. 125 In this work, the electrolyte potentials were also references against the experimental lithium redox voltage, although a slightly different value of 1.37 V was used. A similar proximity effect on electrochemical stability was found by Rajput et al. 129 in a study of Mg-electrolytes. They found that the reduced TFSI (1-ethyl-3-methylimidazolium) anion undergoes significant bond weakening when paired with a Mg cation in the solution.
The decomposition of electrolyte molecules upon reduction at the anode surface is a key process during the formation of the solid-electrolyte interface. Wang et al. 116 investigated the reductive decomposition and polymerisation of ethylene carbonate (EC) with first principles calculations, also employing the PCM solvent and evaluating the redox potentials with the method by Zhang. The calculated reaction mechanisms show that a reduction reaction involving multiple EC molecules (2-4) has a lower reduction potential than the single-molecule reaction, indicating that the accurate calculation of redox potentials may require some explicit solvent molecules. The decomposition products predicted by Wang et al. are in good agreement with experimental observations (see references within ref. 116) including the release of ethylene upon reaction.
Leung studied the reductive decomposition of EC on the surface of a graphite anode with direct AIMD simulations (see also section 'Direct simulation of the diffusion dynamics'). 130 On the basis of an AIMD trajectory over 7 ps, the author identifies two reduction pathways resulting in the release of either ethylene or carbon monoxide, both of which had previously been observed in experiments. In contrast to the work by Wang et al. the AIMD simulations by Leung did not impose the reaction mechanism in the computation.

CHALLENGES AND PERSPECTIVES
The overview provided in the previous sections shows that many key properties of lithium batteries, such as the voltage, rate capability and thermal stability, can be reliably addressed by first-principles calculations. Over the last two decades, the accuracy of ab initio methods (thanks to DFT+U and hybrid functionals) and the robustness of the computational methodologies have matured to the point that simulations can often be completely automated, enabling, for example, high-throughput processing of large structural databases. 66,[131][132][133][134] Notwithstanding these advances, every DFT calculation requires as input a structure model in form of atomic positions, and any computation becomes meaningless if the researcher has no clear conception of the relevant materials phases or makes invalid assumptions about the atomic structure. Today, computational battery research is therefore less limited by technical challenges than by our insufficient understanding of the active phases and relevant processes in some materials.
This challenge has become more pressing since the advent of lithium excess cathode materials which often exhibit (partial) cation disorder and undergo non-coherent phase transitions upon charge and discharge. 113,135 Some lithium excess materials have furthermore been argued to gradually change their stoichiometry due to loss of oxygen gas. 136 Although thermodynamically stable phases can generally be discovered with the methodology of section 'Thermal stability of cathode materials', it is still challenging to predict metastable phases that may be present under operation conditions, such as nanostructures, 137 polymorphs 138 or disordered phases 113 (kMC simulations discussed in section 'Diffusivity based on 0 K migration barriers' are one option). However, without a clear understanding of the structural evolution of a material, it is virtually impossible to make quantitative computational predictions, and it is extremely challenging to identify the material's performance limiting attributes. The difficulty to predict the formation of metastable phases also hampers the computational design of entirely new materials: Presently, we are simply not very good in predicting whether a hypothetical material is actually synthesisable.
Most of the first-principles work on battery materials so far has been focused on crystalline solids, owing to the great interest in crystalline cathode materials but also to the challenges involved in first principles modelling of unordered and molecular materials. However, as touched on in the section 'Electrochemical stability of electrolytes', liquid molecular electrolytes are a crucial element for lithium-ion battery safety and rate capability, and the non-crystalline solid-electrolyte interface has an important role for the battery performance. 7 The simulation of such non-ideal (non-periodic) structures typically requires large length scales that are beyond the means of present first principles methods, but may be accessible by coarse-grained or empirically parametrised models, such as classical force-field-based molecular dynamics simulations 139 or phase-field models. 140 Systematic first-principles computational work in the area of molecular electrolytes is only emerging, but recent initiatives, such as the Electrolyte Genome Project, 141 promise to deliver an increased momentum in the next years.