Introduction

Rechargeable lithium-ion batteries (LIBs) have been a breakthrough solution for many energy storage applications. However, we have reached a fundamental limit of what we can achieve with the existing LIB technology1,2, with further advancements demanding high-energy density, improved safety, and reduced costs, while maintaining a high rate performance. For instance, the current target of the US Advanced Battery Consortium for batteries used in electric vehicles is to achieve an 80% state-of-charge in 15 min3,4. Recent advancements in using lithium metal anodes, along with compatible solid-state electrolytes can achieve energy densities significantly beyond current LIBs, but poor charging rates (coupled with growth of lithium dendrites) hamper their performance5. Multivalent (MV) batteries6,7,8, such as those involving Mg ions, provide an alternative to LIBs in offering systems with higher volumetric energy densities (due to the use of a metallic anode9,10,11), but also suffer from poor rate performance7,12. Thus, it becomes essential to understand and improve rate performance in batteries, where computational tools can play a key role.

A critical factor that influences the rate performance of batteries is the diffusion of the electroactive ions in the active materials (i.e., solid electrodes and/or electrolytes), which in turn is largely influenced by the ionic migration barriers (Em) within host frameworks13. The conventional theoretical methodology to estimate Em is using the computationally intensive density functional theory14,15 (DFT)-based nudged elastic band16 (NEB) calculations and/or ab initio molecular dynamics (AIMD) simulations17,18. Our work focuses on DFT-NEB calculations, which usually is the preferred mode of estimating Em, especially in frameworks with “high” barriers19. This is because the accuracy of AIMD calculations depends on the robust sampling of migration events17, which can be difficult in events with low jump frequencies (caused by high Em). Note that the central approximation in performing either DFT-NEB or AIMD calculations is the choice of the exchange-correlation (XC) functional, which broadly describes the nature of bonding and the quantum mechanical interactions among the electrons present in a system.

A popular choice of XC functional (for DFT-NEB and AIMD) is the Perdew–Burke–Ernzerhof20 (PBE) generalized gradient approximation (GGA), due to its reasonable trade-off between computational cost and accuracy21,22,23,24. However, PBE suffers from several shortcomings, such as self-interaction errors (SIE) in transition-metal (TM) systems25,26, overbinding of the oxygen gas molecule27,28,29, and the inability to account for van der Waals interactions associated with layered frameworks30,31,32. Several modifications have been made to improve upon PBE to capture specific ground state properties of interest, resulting in functionals, such as PBEsol33, PBEfe34, and AM0535, but none have fully replaced PBE. Another approach to mitigate inaccuracies incurred within the general class of GGA functionals (including PBE) is to employ a meta-GGA framework, such as the Minnesota36 or the strongly constrained and appropriately normed (SCAN37,38) functional. Importantly, the SCAN functional satisfies all the 17 known constraints of an XC functional (PBE does not satisfy)39, does not overbind oxygen molecule40, accurately predicts formation enthalpies of main group compounds41, estimates the ground state structural and electronic properties better (than GGA)42, and usually predicts the correct ground state polymorph in compositions where multiple polymorphs can exist43. As SCAN has been developed recently, it has not been widely benchmarked and used to calculate Em, yet, in solids. Thus, one of the aims of this article is to test the accuracy of SCAN in Em predictions in materials for energy storage.

In the case of TM-containing oxide frameworks, both PBE and SCAN exhibit significant SIEs26,40,44,45. One of the techniques used to reduce the SIE, in both PBE and SCAN, is to introduce an empirical Hubbard U parameter46, resulting in a PBE+U44,47 (or equivalently GGA+U) and SCAN+U40 functional, respectively. While the computational cost for a GGA+U or SCAN+U is similar to the corresponding GGA or SCAN calculation, GGA+U/SCAN+U is better at predicting the electronic and redox properties of TM-containing compounds compared to GGA/SCAN40,44,48,49. Nevertheless, the Hubbard U-corrected functionals are not widely used for estimating Em in solids, primarily because of convergence difficulties in performing NEB calculations50,51. Additionally, U-corrected functionals are expected to overestimate Em, particularly in electrode frameworks that form polarons52,53,54 within the bulk structure. This is because the migration of a cation in such structures involves the migration of a localized electron, creating configurations with partial electronic occupancies across different TM sites, especially “near” the transition state (TS). Such partially occupied states are penalized by the Hubbard U framework, resulting in an underestimation of the TS energy and an overestimation of Em. Nevertheless, a robust evaluation of the computational performance and accuracy of GGA and SCAN and their corresponding U-corrected versions is required before any specific computational framework can be employed in a more widespread manner.

To date, studies have focused on understanding the influence of XC frameworks, especially with SCAN, on several “static” properties, such as intercalation voltages55,56,57,58, band gaps55, interlayer spacings55, and magnetic properties57, and not on “dynamic” properties such as Em. In this study, we perform a detailed investigation of the impact of the XC functional choice on Em predictions in six different well-studied electrodes, namely the layered-LiCoO2 (space group: \(R\overline{3}mH\)), the spinels-LiMn2O4 (\(Fd\overline{3}m\)), MgMn2O4 (I41/amd), and MgxTi2S4 (\(Fd\overline{3}m\)), the olivine-LiFePO4 (Pnma), and the post-spinel-NaV2O4 (Pnma). Also, we consider three different solid electrolytes, namely the spinel-MgSc2Se4 (\(Fd\overline{3}m\)), the tetragonal-Na3PS4 (\(P\overline{4}{2}_{1}c\)), and the orthorhombic-Li3PO4 (Pnma). The choice of materials studied in this work are motivated by: (i) the availability of experimental mobility data, (ii) the diversity of structural frameworks (layered, spinel, post-spinel, etc.) with different redox-active TMs (such as, Co, Mn, Fe, V, Mn, and Ti), and (iii) the heterogeneity of intercalant (or mobile) cations, including Li+, Na+, and Mg2+.

We consider four different XC functionals for calculating Em in electrode materials: GGA (PBE), SCAN, GGA+U, and SCAN+U. While a number of strategies have been proposed in determining U values59, we employ U corrections, with both GGA+U and SCAN+U, that are optimized to reduce errors against experimental oxidation enthalpies40,44,45,60. This choice is ideal since ionic migration in electrodes often involves a “local” redox process, with the migrating ion typically causing the movement of electron(s) as well50,61. Furthermore, modeling ionic migration via vacancy-mediated-hops results in the presence of mixed TM oxidation states, which are best described with an “averaged” U parameter. Apart from the XC functionals considered, we have also analyze the role of adding a uniform background charge (UBC) and using the climbing image (CI16) approximation on calculated Em in the solid electrolytes considered (see Migration barriers section under Results).

The main objective of our work is to validate the calculated Em with the various XC approximations with existing experimental data for the nine prototype materials mentioned above. Overall, we observe that SCAN provides a marginally better quantitative accuracy of Em across all materials considered on average, compared to GGA and the U-corrected frameworks, albeit with several exceptions. In solid electrolytes, we do not find the addition of a UBC and/or the CI framework to significantly alter Em with both GGA and SCAN, except for Li3PO4. Importantly, Em evaluated with SCAN are larger than those with GGA across several materials, which leads to SCAN (GGA) providing a qualitative upper (lower) bound to the experimental Em. Furthermore, the addition of U corrections (in electrodes) has a lower impact on SCAN-Em than GGA, which can be attributed to the better description of the underlying electronic structure, as characterized by the on-site magnetic moments on TMs. Therefore, among the XC functionals considered here, we expect SCAN to be the best choice of XC framework for Em predictions, although with higher computational costs and convergence difficulties than GGA-based functionals. Also, we note that GGA can provide better Em estimates than SCAN in specific materials, such as MgxTi2S4, MgSc2Se4, and Li3PO4, thus making Em evaluations with GGA quite useful. Our work provides a better understanding of the underlying physics behind the estimation of Ems in energy storage materials, and will aid in the selection of the right XC functional for a given structure in future studies, thus enabling the discovery of ion conducting electrodes and solid electrolytes via computational workflows62.

Results

Structures

Figure 1 shows the structures of the materials considered in this work, with additional details reported in Supplementary Table 1 of the electronic supporting information (ESI). Li3PO4, which is a solid electrolyte with high Li-conductivity and ease of preparation63 can be experimentally synthesized in two different phases, γ and β64. The γ-Li3PO4 polymorph (Pnma space group, Fig. 1a) exhibits oxygen atoms arranged in a hexagonal-close-packed sub-lattice, with Li and P atoms occupying 16 out of 32 available tetrahedral voids. We consider the γ polymorph instead of β since it is an ordered configuration of Li3PO4 despite being thermodynamically metastable65. Similar to Li3PO4, the Na ionic conductor, Na3PS4 also exhibits one of two polymorphs, the tetragonal-Na3PS4 (or t-Na3PS4, \(P\overline{4}{2}_{1}c\), (Fig. 1b)) and the cubic-Na3PS4 (or c-Na3PS4, \(I\overline{4}3m\)), at 298 K, depending on the synthesis technique66,67,68. The main difference between the two structures of Na3PS4 is the split of the Na(6b) site in the cubic phase into two sites, Na1(2a) and Na2(4d) in the tetragonal phase67,68, as indicated in Fig. 1b. We investigate the t-Na3PS4 as it is stable and exhibits an ordered structure69.

Fig. 1: Supercells of all the structures considered in this work.
figure 1

The migrating and the transition-metal ions are represented by orange and green spheres/polyhedra respectively. Phosphate ions are shown by purple polyhedra.

Full size image

Figure 1d displays the O3 phase70, an important ordering of LixCoO2, which exhibits oxygens in a “ABCABC”-type stacked close-packed layers with the Li and Co positioned in octahedral sites across alternate cation layers71. We consider the so-called divacancy (DV72) hop for estimating the Em in LixCoO2, while the calculated Em for the monovacancy hop is included in Supplementary Fig. 1 of the ESI. Another commercial electrode material for LIBs is LiFePO4 (Fig. 1h) that exhibits an ordered olivine-type structure (Pnma), consisting of corner-sharing FeO6 octahedra and PO4 tetrahedra with the Li ions residing in “tunnels" along the b-axis. Both LiCoO2 and LiFePO4 are well-studied Li-ion cathode frameworks73,74. The Na-cathode considered in this work is the post-spinel NaV2O4 (Pnma, Fig. 1e), where the mobile Na ions are located in the center of 1D tunnels formed by the edge sharing VO6 octahedra, along the b axis.

Among the four spinels considered here, the MgMn2O4 (Fig. 1g, I41/amd) is one of the Mg oxide-cathodes to exhibit a reasonable extent of Mg mobility75,76. MgMn2O4 adopts either a tetragonal or cubic symmetry depending on the concentration of Jahn-Teller distorted Mn3+, which is determined by both the Mg content and the degree of spinel “inversion”75. Ideally, for xMg = 1, the structure of MgMn2O4 is tetragonal. On the other hand, spinels-MgSc2Se4 and LiMn2O4 (Fig. 1f) are cubic (\(Fd\overline{3}m\)) since they have cations that do not exhibit cooperative Jahn-Teller distortion. In the case of spinel-MgxTi2S4 (Fig. 1c), the intercalated Mg atoms occupy octahedral sites within the cubic structure instead of the typically occupied tetrahedral positions77. Note that we have calculated the Em for MgxTi2S4 in the high-vacancy limit (i.e., x → 0), since the Mg sub-lattice is disordered at high Mg concentrations (at x → 1). For all electrodes except MgxTi2S4, we have calculated the Em at the dilute-vacancy limit (see Methods section).

Migration barriers

We compare barriers obtained from the various XC frameworks (colored bars) against experimental data76,77,78,79,80,81,82,83,84 (horizontal dashed lines) of electrodes and solid electrolytes in the top and bottom panels of Fig. 2. The detail on synthesis and measurments of ionic mobility employed experimentally for all materials are listed in Supplementary Table 2 of the ESI. Numerical values of Em are listed in Table 1. Minimum energy pathways (MEPs) and Em not included in the main text are compiled in the ESI (Supplementary Figs. 1–8). In solid electrolytes, we do not consider the U-corrected functionals since these materials do not contain redox-active TMs with partially filled d orbitals85. Note that ionic migration in solid electrolytes considered here occurs via the vacancy-mediated mechanism, which requires the creation of a Li/Na/Mg vacancy within the stoichiometric supercell used in the NEB calculation. Subsequently, a UBC is often added to the periodic cell (referred to in this work as “NELECT”)86,87 to maintain charge neutrality during the calculation. Since the impact of adding a UBC on Em has not been rigorously quantified in diverse solid electrolytes before, we have calculated Em using both GGA and SCAN XC with and without NELECT. Note that interstitial/interstitialcy-based migration mechanisms can be active in solid electrolytes, particularly in disordered systems88, but not considered in this work.

Fig. 2: Computed barriers, using different XC frameworks, for the electrodes (top panel) and solid electrolytes (bottom panel) considered in this work.
figure 2

Horizontal dashed lines indicate experimentally reported values. “ne” in the electrolytes panel indicates the inclusion of UBC (i.e., NELECT) in the NEB calculation, while bars corresponding to “CI” highlight the usage of the climbing image framework. Barriers plotted for Li3PO4 are for the 1(d) → 2(c) hop. “DV” in LiCoO2-DV stands for the divacancy hop mechanism.

Full size image
Table 1 Experimental and computed migration barriers Em (in meV) of all materials considered in this work.
Full size table

The CI framework involves a modification to the regular NEB in removing the spring forces surrounding the image with the highest energy, nominally allowing for a better description of the TS (i.e., the saddle point)16. While CI-NEB has been shown to more accurately predict Em, it can be computationally difficult to converge89. Hence, for NEB calculations in solid electrolytes, we also perform a CI-NEB calculation for each XC framework considered, after the regular NEB converged, to estimate any changes in the calculated Em.

An examination of the data suggests that the magnitude of Em predicted by SCAN is higher compared to GGA, while the barriers reported by SCAN+U are lower than GGA+U in electrodes (with NaV2O4 being the only exception). In solid electrolytes, the addition of NELECT does not impact the GGA or the SCAN barriers as seen in MgSc2Se4 and Na3PS4. In contrast, adding NELECT does affect the barriers of Li3PO4, by ~9% for GGA and ~4% for SCAN (for the 1(d) → 2(c) hop, see Li3PO4 section for additional details). In general, adding NELECT reduces Em marginally. The addition of CI in our calculations does not change the barrier for a given functional (except in the case of Li3PO4), with variations of Em < 2% across all corresponding XC functionals, with and without CI.

Overall, the barriers predicted by SCAN have a better numerical accuracy compared to GGA or the Hubbard U-corrected frameworks, with a mean absolute error (MAE) of ~140 meV when compared to 178.33, 183, and 145.17 meV for GGA, GGA+U, and SCAN+U, respectively. There are specific exceptions where GGA performs better than SCAN including NaV2O4, MgxTi2S4, MgSc2Se4, and Li3PO4. Similarly, in LiFePO4, LiMn2O4, and NaV2O4, GGA+U exhibits lower errors (MAE ~177 eV) than SCAN+U (~212 eV). Further, we discuss in detail systems, including γ-Li3PO4, t-Na3PS4, spinel-MgMn2O4, and layered-LiCoO2, in the following sections to highlight the general trends and specific anomalies.

γ-Li3PO4

Out of the several possible migration paths along the three orthogonal axes of Li3PO4, the migration of a Li-vacancy occurs along a zig-zag path across the two nonequivalent Li sites, namely Li1 and Li2, which are also referred to as “d” and “c” sites, as in Fig. 1a. Specifically, the presence of these two Li sites leads to a 1(d) → 2(c) → 3(d) hop, requiring the presence of two metastable Li vacancies, as shown by the blue spheres in Fig. 3a and the black arrows in its inset. Note that 1(d) and 3(d) are symmetrically equivalent Li1 sites, but are named so to distinguish them during Li migration, as per the notation by Du and Holzwarth90. The energy of formation of a Li vacancy over the d and c sites, with multiplicity 8 and 4, respectively, differ by ~0.22 eV90. There are three Li-hops that can be calculated in Li3PO4: (i) 1(d) → 2(c), (ii) 2(c) → 3(d), and (iii) 1(d) → 2(c) → 3(d).

Fig. 3: Migration pathway and computed barriers for γ-Li3PO4.
figure 3

a 1 × 2 × 2 supercell of Li3PO4 used for the NEB calculations. Green and purple polyhedra represent LiO6 and PO4 groups, respectively. The inset shows the migration of the Li ion along the 1(d) → 2(c) → 3(d) path (blue spheres). b MEPs and Em are plotted for the 1(d) → 2(c) pathway, without NELECT. GGA-CI and SCAN-CI indicate CI-NEB calculations with GGA and SCAN, respectively. c MEPs and Em for the 1(d) → 2(c) path, with NELECT added.

Full size image

Figure 3b, c show the MEPs and the corresponding barriers (in meV) for the 1(d) → 2(c) hop, without and with NELECT, respectively. The barriers pertaining to the 1(d) → 2(c) → 3(d) divacancy and 2(c) → 3(d) hops are included in the ESI (Supplementary Fig. 7). In general, Em predicted by SCAN XC functional in Li3PO4 (for all hops considered) are higher than the corresponding GGA approximation. The addition of CI does not alter the magnitude of Em significantly, except in the case of the 1(d) → 2(c) → 3(d) hop which shows a maximum deviation of ~16% (between GGA and GGA + CI, without NELECT). The qualitative nature of the MEPs captured by both GGA and SCAN XC functionals is similar, and does not change significantly with the inclusion of NELECT. However, adding NELECT does change the magnitude of Em significantly in some hops of Li3PO4, with the variation higher in GGA-based frameworks (~3–35%variation from the corresponding Em without NELECT) than SCAN-based (~2–23% variation).

The barrier for the 1(d) → 2(c) hop is lower than for the 2(c) → 3(d) hop, which is consistent with the report by Du and Holzwarth90 with our calculated GGA-CI and GGA Em also exhibiting good agreement (see Supplementary Fig. 7). The Em for the 1(d) → 2(c) → 3(d) hop, when considered along with the vacancy formation energy is consistent with the effective activation barrier reported in the impedance spectroscopic study by Ivanov-Shitz et al.91. The experimentally reported Em for Li+ migration in a solid solution of Li3+xP1−xSixO4 (0 < x < 0.4) with the γ-Li3PO4 structure, via nuclear magnetic resonance (NMR) measurements, is ~27.6 ± 0.7 kJ per mol (~286.05 ± 6 meV) for 80 mol% Li3PO484. The Em predicted by GGA+NELECT for the 1(d) → 2(c) hop is closest to the experimental value (albeit a significant ~31% error) while SCAN+NELECT overestimates it by ~ 52%. The addition of NELECT does improve the numerical accuracy of both GGA and SCAN (for the 1(d) → 2(c) hop) with respect to experimental data. Note that our on-site magnetic moment data does not indicate any significant electron localization upon NELECT addition (with either GGA or SCAN), ruling out any localization of electrons contributing to better Em agreement with experiments.

Tetragonal-Na3PS4

The t-Na3PS4 supercell used in our NEB calculation is shown in Fig. 4a, with the blue spheres indicating the migrating Na atom. Orange and yellow spheres represent the Na1 and Na2 sites, respectively. Purple squares represent the PS4 polyhedra. There are two pathways in which the Na ion can diffuse through the structure: (i) hop from Na2 to Na1 site along with the a or b axis and (ii) hop from Na2 to a neighboring Na2 site along the c axis. We investigate the Na1 → Na2 hop (Fig. 4a) since a previous study has reported that the Na2 → Na1 hop is more facile than the Na2 → Na2 hop83. Also, several experimental studies have been conducted to study the Na ion mobility in Na3PS4. For example, Chu et al.92 reported that pristine t-Na3PS4 exhibits a conductivity of ~0.01 mS/cm at 300 K, corresponding to an Em of ~317 meV, using impedance spectroscopy, while similar measurements by Krauskopf et al.93 and Chang et al.83 reported Em of ~380 and ~350 meV, respectively.

Fig. 4: MEPs and the corresponding barriers for Na3PS4.
figure 4

a 2 × 2 × 2 supercell of Na3PS4. Na1 and Na2 sites are represented by orange and yellow spheres, respectively. Blue spheres indicate the Na migration pathway. b MEPs and Em for the Na1 → Na2 hop. GGA + NE and SCAN + NE indicate calculations with NELECT added.

Full size image

Figure 4b shows the calculated MEPs and Em for the GGA and SCAN functionals, with and without NELECT. The magnitude of Em is unaltered on including the CI approximation and are only marginally affected by the inclusion of NELECT with both GGA and SCAN (maximum change of ~8 meV). The shape of the MEPs are qualitatively similar across all the four XC frameworks. Similar to Li3PO4, SCAN-based calculations predict a higher Em compared to corresponding GGA-based versions. Moreover, the difference in Em with NELECT addition in SCAN (~3 meV) is lower than in GGA (~8 meV), similar to trends observed in Li3PO4. However, all our calculated Em are significantly lower than experimental reports (by ~300 meV), while classical MD simulations have predicted barriers in better agreement with experiments 94.

Note that our calculated Em are consistent with the GGA-Em of Huang et al.66 and Yu et al.95 for the Na2 → Na1 hop, suggesting that the large discrepancy between DFT and experimental activation energy values may be largely attributed to the non-inclusion of defect formation energies in Em calculated by NEBs. Indeed, previous AIMD simulations66,96 have indicated a decrease in Em (from >300 meV in pristine-Na3PS4 to <280 meV) due to the presence halogen dopants (which nominally create Na-vacancies) in both c-Na3PS4 and t-Na3PS4.

Spinel-MgMn2O4

Spinel-MgMn2O4 has its oxygen anions adopting a face-centered-cubic lattice, with half of the octahedral voids occupied by Mn and 1/8th of the tetrahedral voids occupied by Mg. Mg2+ in MgMn2O4 migrate across a 8a → 16c → 8a pathway, where the initial and final 8a sites share a triangular face with the intermediate, empty octahedral (16c) site. The edges of the 16c site are shared by six other Mn octahedra, resulting in a “ring” of Mn sites adjacent to it75. Figure 5a provides a closer look of the 16c and the adjacent ring sites (solid yellow lines), with the left and right images indicating the view from “top” and “bottom”, respectively. The text annotations on each Mn site in Fig. 5a identify the specific Mn ring site. Green and purple octahedra, respectively, represent the initial and the final 8a sites of the migrating Mg ion. The yellow upright (inverted) triangle signifies Mn octahedra coming out (pointing into) the plane of the paper. The inverted triangle represents Mn24 and Mn23 in the left and right subpanel, respectively, while the upright triangle signifies Mn20 in both subpanels. Another view of the Mn sites adjacent to the initial and final 8a sites is provided in Supplementary Fig. 9a.

Fig. 5: Migration barriers computed for spinel-MgMn2O4.
figure 5

a Top and bottom views of the Mn “ring” sites (dashed yellow lines) adjacent to the intermediate 16c site during Mg migration. The initial and the final Mg ions are represented by the green and purple tetrahedra, respectively. The yellow triangles represent additional Mn sites (see text for further description). b MEPs and Em. c On-site magnetic moments of Mn ions in the MgMn2O4 structure as a function of image numbers (i.e., as Mg is undergoing migration), where 1st and 9th images correspond to the initial and final Mg sites, respectively. The number indicated on the TM sites in (c) is identical to the number indicated in text labels for Mn sites in (b).

Full size image

Panel b of Fig. 5 shows the MEPs (left subpanel) and the corresponding Em as bar charts (right subpanel). Figure 5c shows the on-site magnetic moments of all TM ions (i.e., Mn) in the MgMn2O4 structure used for the NEB calculation, as a function of the image number (where each image represents a snapshot of Mg during its migration). Specifically, image number 1 (9) represents the initial (final) Mg site. The TM number in panel c is identical to the notation used for labeling Mn ions in panel a. The on-site magnetic moments are represented via the heat map ranging from 2.5 μB to 4 μB, where the lower (higher) values represent Mn4+ (Mn3+). Importantly, changes in magnetic moments track the movement of electrons as Mg migrates within the spinel.

Note that due to convergence difficulties, we had to use a different supercell and relaxation scheme for the GGA-NEB calculation of MgMn2O4, which explains the different number of TMs for GGA (see Supplementary Fig. 9b). Notably, GGA wrongly predicts the bulk cubic version of spinel-MgMn2O4 to be the ground state instead of the tetragonal, Jahn-Teller distorted structure97. This is because GGA wrongly delocalizes the valence electrons between Mn and O atoms to an extent that it cannot account appropriately for the Jahn-Teller distortion of Mn3+. In contrast, SCAN(+U) and GGA+U functionals, with their intrinsic tendency to decrease the SIE, rightly predict the tetragonal polymorph to be the ground state. Therefore, care must be taken to ensure that the “right” structure is used for calculating Em with GGA in the case of Jahn-Teller distorted structures. To address this issue here, we have used a PBEsol+U-relaxed spinel conventional cell containing 56 atoms (without Mg vacancy) to initialize the endpoint images and subsequently the NEB calculation, which is identical to the procedure of Gautam et al. 75.

From Fig. 5b, we observe that the MEPs are qualitatively similar, with all functionals displaying a local energy minimum corresponding to the 16c site. Importantly, the Em calculated with GGA+U, SCAN, and SCAN+U tend to be in reasonable agreement with the experimentally reported Em (700 ± 60 meV)76. Specifically, the maximum and minimum error, with respect to experiment, among these XC frameworks is about 13% (SCAN) and 4% (SCAN+U), respectively. On one hand, the magnitude of error reported experimentally (i.e., a window of 120 meV) is ~17% with respect to the 700 meV average value, which is approximately the same order of magnitude as the error among the predictions of GGA+U, SCAN, and SCAN+U, highlighting the normal magnitude of errors that can be expected across theory and experiments for electrodes. On the other hand, GGA does underestimate the experimental Em significantly by ~32%. Thus, the addition of the U correction increases the numerical accuracy significantly in the case of GGA, while it does not make a drastic difference in the case of SCAN for MgMn2O4. The smaller difference between SCAN and SCAN+U Em predictions, as compared to GGA and GGA+U, can be due to the better description of XC interactions and lower SIEs in SCAN.

Similar to other materials considered in this work, SCAN does predict a higher Em than GGA. However, SCAN+U’s barrier is lower than GGA+U’s, which may be attributed to the higher magnitude of U correction used in GGA+U (3.9 eV60) as compared to SCAN+U (2.7 eV40). All four XC functionals predict changes of on-site magnetic moments on Mn atoms as the Mg2+ migrates, indicating a concurrent migration of electrons resulting in a polaronic hop that is characteristic of semiconducting/insulating electrodes53. The expected behavior is that as the Mg2+ ion hops from the initial to the final 8a site (i.e., as image number increases in Fig. 5c), the oxidation state of the Mn ions closer to the initial 8a site will increase from +3 to +4 (corresponding to a decrease in their magnetic moments), and the vice-versa for Mn ions that are closer to the final 8a site.

Importantly, we find that changes of magnetic moments (and hence electron hopping) are more accurately captured by SCAN, SCAN+U, and GGA+U than plain GGA, as indicated by two distinct Mn atoms displaying lower magnetic moments (i.e., Mn4+ oxidation state), corresponding to the Mg-vacancy migrating, throughout the NEB. For instance, in SCAN+U, the holes (Mn4+) are located on Mn12 and Mn14 sites that are closer to the final Mg site (purple tetrahedron in Fig. 5a) at the beginning of Mg migration (i.e., holes are closer to the vacancy). Subsequently, as Mg migrates to its final position, the holes also migrate to Mn4 and Mn6 sites that are closer to the initial Mg site (green tetrahedron). Thus, SCAN+U captures the hole (electron) hop that occurs alongside the Mg hop. In the case of SCAN, Mn20 tends to have its magnetic moment constant throughout the NEB as the Mn20 octahedron shares its vertex with both the initial and the final positions of the Mg tetrahedron, and hence accommodates the hole (Mn4+) better as the Mg-vacancy migrates. Mn24, which is closer to the initial 8a, has its magnetic moment decreased as Mg moves from image 1 to 9 in the SCAN-NEB calculation, while the opposite behavior is observed in the case of Mn23, which is closer to the final 8a.

Layered-LiCoO2

The diffusion of Li ions in layered-LiCoO2 during the charging and discharging process occurs across each Li layer. The two possible Li migration pathways are: (i) the tetrahedral site hop (TSH) or divacancy hop, and (ii) the oxygen dumbbell hop (ODH) or monovacancy hop, as illustrated in panels a and b of Fig. 6, respectively72. Hollow black rectangles in Fig. 6 represent Li vacancies in the octahedrally coordinated Li layers. The arrows correspond to the migration path of the Li-ion, specifically, via the intermediate tetrahedral site in the TSH mechanism and passing through an O–O bond in ODH. Previous studies have shown that Li ions prefer the TSH mechanism, involving the presence of an additional Li-vacancy adjacent to the migration pathway (hence termed DV hop) instead of the ODH mechanism, which can be attributed to electrostatic repulsions of the migrating Li with other Li and Co cations71,72.

Fig. 6: MEPs and the corresponding barriers for LiCoO2.
figure 6

a, b display the divacancy and monovacancy Li migration mechanisms, respectively, in LiCoO2. c Calculated MEPs and Em. d On-site magnetic moment changes of the Co ions as a function of Li-ion migration with the GGA+U (top panel) and GGA/SCAN/SCAN+U functionals (bottom panel).

Full size image

While Fig. 6c plots the MEPs and Em associated with the DV hop, the data for ODH is provided in Supplementary Fig. 1. The MEP of the DV hop shows a local minima corresponding to the intermediate tetrahedral site for all functionals, and the calculated Em are lower than the monovacancy hop, as expected. The experimental barrier for LiCoO2 is ~300 meV, as observed in NMR studies of Nakamura et al.78. Notably, SCAN predicts a barrier that is identical to the experimental value, while GGA underestimates the barrier by ~33%. Both GGA+U and SCAN+U overestimate Em versus experiment, by ~73% and 26%, respectively. However, the qualitative trends remain identical across all the XC functionals. Note that due to convergence difficulties, we performed the SCAN+U-NEB calculation using relaxed SCAN-NEB structures (see Section S1 of ESI for details).

Figure 6d displays the on-site magnetic moments of the TM (i.e., Co) ions in the LiCoO2 supercell as a function of the image number. The magnitude of the heat map ranges from 0 to 3 μB, where lower (higher) values indicate Co3+ (Co4+). Unlike MgMn2O4, only GGA+U tends to capture a local electron hop, as indicated by changes in magnetic moments across several Co ions as Li migrates, whereas the other functionals exhibit negligible changes in magnetic moments. Note that LixCoO2 is a metallic compound over a large range of Li concentrations71, which indicates that the magnetic moments of Co ions should not change significantly as Li migrates (due to electron delocalization in a metal). This aspect is captured well by all functionals except GGA+U, even though all four functionals predict semiconducting behavior in pristine-LiCoO2, without any Li-vacancy45. Thus, the large error in the GGA+U-calculated Em (versus experiment) can be attributed to the spurious localization of electrons by GGA+U in a metallic structure.

Discussion

Given that good mobility of the intercalating cation in electrode and solid electrolyte structures and a proper computational framework for estimating ionic mobility with reasonable accuracy is essential for developing high-energy-density battery systems, we analyze the accuracy of a few commonly employed XC approximations within DFT in predicting Em in battery materials. Specifically, we benchmark predictions of GGA, SCAN, GGA+U, and SCAN+U XC frameworks for six distinct electrodes, and consider GGA and SCAN (with and without NELECT/CI) for three topical solid electrolytes. Importantly, we observe that SCAN has better numerical accuracy in estimating Em across the structures considered than other XC frameworks. Additionally, we observe that the difference in Em predicted by SCAN and SCAN+U (MAE across electrodes is ~51 meV) is consistently lower than the difference in Em predicted by GGA and GGA+U (MAE ~200.5 meV). Thus, for Em predictions, we do not expect the addition of U to make a significant difference with SCAN. This behavior can be attributed to overall lower SIEs in SCAN than GGA. The addition of NELECT or CI does not affect the qualitative or the quantitative nature of the barrier in the case of GGA or SCAN in solid electrolytes (except in Li3PO4).

In electrodes, GGA+U, SCAN, and SCAN+U XC frameworks capture changes of on-site magnetic moments, which quantify the extent of electron hopping that occurs in conjunction with an ionic hop, better compared to GGA, in the case of semiconducting MgMn2O4. However, in the case of metallic LixCoO2, we find GGA, SCAN, and SCAN+U to be better than GGA+U. Thus, the specific choice of the XC functional, especially between GGA and GGA+U, will depend on the electrode’s electronic nature, while we expect the addition of U (or not) to make a less significant impact on Em calculations with SCAN, irrespective of the underlying electronic structure.

Given that U is a “tunable” parameter in GGA+U and SCAN+U calculations, it is important to examine the effect of changing the U magnitude in both XC functionals in predicting Em. We have calculated the Em for LiMn2O4 with GGA+U and SCAN+U, where the U values are 2.7 and 3.9 eV, respectively, with the results plotted in Supplementary Fig. 5 of the ESI. Note that the U values optimized for oxidation enthalpies for Mn-oxides are 3.9 and 2.7 eV with GGA+U and SCAN+U, respectively40,60. Thus, we vary the U magnitude by ~1 eV, which is usually sufficient to observe qualitative differences in predictions. Furthermore, we observe that the U value, derived by from linear-response theory47, for MnO (averaged over its rocksalt and zinc blende polymorphs) with SCAN+U is 2.8 eV98, which is similar to the optimized U used in this work. For GGA+U, the value derived with linear-response (3.2 eV) is within the range of U values (2.7–3.9 eV) explored for LiMn2O4 in this study98.

On one hand, we find that changing the U value does not change the Em predictions with GGA+U (349-364 meV, see Supplementary Figs. 4, 5), with values at U = 2.7 and 3.9 eV in reasonable agreement with experiment (~500 meV, Table 1)79. On the other hand, SCAN+U predictions of Em vary significantly with varying U (327–750 meV), with the SCAN+U = 2.7 Em (327 meV) in better agreement with experiment. The GGA+U = 3.9 and SCAN+U = 2.7 NEB calculations capture better changes of the on-site magnetic moment of the Mn sites as the Li+ migrates in LiMn2O4 (see Supplementary Fig. 5). Thus, we can conclude that SCAN+U’s Em predictions are more sensitive to changes in U values compared to GGA+U. We expect the U values optimized for oxidation enthalpies to describe better the underlying electronic structure with both XC approximations.

Apart from SCAN-calculated Em being higher than GGA-calculated Em in general, SCAN overestimates the experimental Em in select materials, such as MgSc2Se4, while GGA underestimates the experimental value in such materials. Thus, we can consider GGA and SCAN to provide qualitative upper and lower bounds to the experimental Em, which will be an useful exercise to quantify uncertainty in materials where experimental measurements of Em can be challenging (such as, decoupling grain boundary contributions or other coupled interactions in impedance spectroscopy measurements). To improve such uncertainty quantifications within computational frameworks, Em calculations using hybrid functionals, especially in electrodes99 can also be attempted. Indeed, Barnes et al.99 demonstrated that a full NEB optimization using hybrid functionals100,101 classified the active and inactive Mg migration pathways in α-MoO3 better than PBE, albeit at a significantly higher computational cost. Another theoretical framework that may yield accurate Em predictions is DFT + U + V, with an additionally tunable “V” parameter that represents inter-site electronic interactions102. Given that experimental errors can be in the ~10% range, any theoretical calculation in a similar ballpark (10–15% error range) can be considered acceptable.

We observe that SCAN has a relatively lower MAE (versus experiments) compared to other XC functionals in Em predictions. However, to decide if SCAN-based NEBs are worth pursuing, considering the marginal gain in MAE, an examination of the computational cost associated with SCAN calculations is necessary. Figure 7 shows the computational time (in minutes) per image in an NEB calculation normalized to the number of cores used per image, and relative to the time taken in a GGA-NEB calculation (which sets the zero on the y-axis). Note that the objective of our computational time comparison is not to rigorously benchmark the performance of DFT codes, but to get a qualitative idea of the relative cost associated with the different XC functionals considered.

Fig. 7: Relative computing time, in units of minutes taken per number of cores used per NEB image for all the Em calculations performed in this work.
figure 7

NE and CI represent addition of NELECT and CI, respectively. Li3PO4(I) and Li3PO4(II) represent 2(c) → 3(d) and 1(d) → 2(c) hop respectively.

Full size image

In general, we find that the computational cost for SCAN (red bars) and SCAN+U (purple bars) are higher than both GGA (blue) and GGA+U (orange) in electrodes, while SCAN appears consistently more expensive than GGA in solid electrolytes (except the 1(d) → 2(c) hop in Li3PO4). However, SCAN does exhibit “similar” computational time as GGA (or GGA+U) in select materials, including, LiMn2O4, MgxTi2S4, MgSc2Se4, and Li3PO4 (2(c) → 3(d) hop). Nevertheless, both SCAN and SCAN+U suffer from convergence difficulties, in addition to increased computational time, especially in electrodes, namely LiCoO2, LiFePO4, and NaV2O4 (see Section S1 of ESI). Shifting to other meta-GGA frameworks, such as the recently developed r2SCAN103, (which improves convergence while retaining SCAN’s numerical accuracy) or the Minnesota functionals36 can be one possible solution to mitigate convergence issues, but rigorous testing is needed. Given the higher computational costs and convergence difficulties of SCAN, GGA may yet be relevant for obtaining “quick” and reasonably accurate Em predictions. Indeed, GGA may even be more accurate than SCAN in some materials, as illustrated by MgxTi2S4, MgSc2Se4, and Li3PO4 in our work.

In conclusion, we have systematically investigated the accuracy of XC frameworks, namely, GGA, GGA+U, SCAN, and SCAN+U in the prediction of migration barriers against existing experimental data for six different electrodes and three different solid electrolytes, spanning diversified crystal systems, to select the appropriate XC candidate for a given structure. Our results indicated that SCAN has better numerical accuracy, on average, than GGA and other Hubbard U-corrected XC frameworks, but with higher computational costs and convergence difficulties. The addition of a uniform background charge or the implementation of the climbing image, in the case of solid electrolytes did not significantly affect the Em predictions of both GGA and SCAN (except for Li3PO4). Also, we found the addition of U corrections to have a lower impact on Em predictions with SCAN than GGA, which can be attributed to a better description of the electronic exchange contribution by SCAN. However, we also observed that changing U values within the SCAN+U framework led to more dramatic changes in calculated Em in LiMn2O4, highlighting the sensitivity of Em on the U value used in the SCAN+U description. Thus, we find that all considered functionals are equally reliable (or unreliable) for a systematically accurate prediction of Em. However, on average, the SCAN functional turns out to provide the most accurate results among the considered functionals, albeit with several exceptions where GGA-based Em predictions are more accurate than SCAN. Note that if qualitative trends in Em are essential and the tolerance towards numerical accuracy (versus experiments) is low, GGA remains a highly useful option for all materials. In cases where experimental Em measurements can be difficult, we observed that GGA- and SCAN-predicted Em can be used as qualitative lower and upper bounds for the experimental Em, since predictions of SCAN-Em are typically larger than GGA. Given that ionic mobility predictions are critical in identifying candidate electrode and solid electrolyte materials for developing high-energy-density battery systems with robust rate performance, our work will help choose the appropriate XC framework for predicting Em in materials.

Methods

Structure relaxations

All DFT calculations were performed using the Vienna ab inito simulation package (VASP 6.1.2) that employs the projector augmented wave (PAW) approximation to describe the effects of core electrons in the many-body wavefunctions104,105,106. The plane wave basis kinetic energy cutoff was set to 520 eV. A fully automatic, Γ-centered grid generation scheme was used to sample the 1st Brillouin zone, at a density of at least 48 k-points per Å for relaxing the bulk structures. For NEB calculations, we reduced the k-mesh density 32 per Å. The input structures for all systems were obtained from the inorganic crystal structure database (ICSD)107. The lattice vectors, cell volume, and ionic positions of all input bulk structures were relaxed until the total energy and atomic forces were below 0.01 meV and ∣0.03∣ eV per Å, respectively. The U values for GGA+U calculations used were identical to the ones used in Materials Project60. We used the U parameters by Gautam et al.40,45 for SCAN+U calculations.

NEB calculations

To estimate Em, the fully relaxed endpoint geometries were obtained, and the MEP was initialized by linearly interpolating both the atomic positions and lattice vectors to create seven images between the endpoints, with a spring force constant of 5 eV per Å between images. The limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) method108 was used to converge to the MEP. The intermediate images were optimized along the reaction coordinate until the total energy converged to within 0.01 meV for each image. The NEB was considered converged if the force component perpendicular to the elastic band was below ∣0.05∣ eV per Å. The barriers were calculated in the dilute-vacancy limit for all the systems, i.e., one Li/Na/Mg vacancy was created in the corresponding supercell with stoichiometric composition51,61,75,109, with the exception of MgxTi2S4. Appropriate supercells of the corresponding primitive cells were used to avoid spurious interaction between the diffusing species by ensuring that the moving ion was always at a distance >8 Å from its periodic image85. All Em calculated in this work correspond to vacancy-mediated mechanisms and are not interstitial-based.

All calculated migration energy profiles are provided either as part of the main text or in the ESI. The structures for which we experienced severe convergence difficulties in our NEB calculations and the strategies we used to overcome such problems are also detailed in the ESI. For Em calculations with UBC in solid electrolytes, we added the number of valence electrons corresponding to the migrating vacancy to the remaining number of valence electrons contributed by the rest of the atoms in the supercell. In VASP the number of electrons is controlled by the keyword “NELECT”, which initializes the total valence charge as a “jellium”-background. For example, our NELECT keyword for a migrating Li vacancy in a Li3PO4 supercell was Σ(valence electrons of remaining Li, P, and O in supercell) + 1. The number of valence electrons for each atom was determined by the corresponding PAW potential employed in the NEB calculation.