Enabling accurate modelling of materials for a solid electrolyte interphase in lithium-ion batteries using effective machine learning interatomic potentials

Abstract

Accurate modelling of the structural and dynamic properties of the solid electrolyte interphase (SEI) in lithium-ion batteries remains a longstanding challenge due to the high complexity of the SEI structure and the lack of structural information. Atomistic simulations using molecular dynamics (MD) can provide insights into the structure of the SEI but require large models and accurate interatomic potentials; however, existing computational tools struggle to evaluate these potentials in mixed-material systems efficiently and reliably. Here, we demonstrate the effectiveness of machine learning interatomic potentials (MLIPs) defined using amorphous structures as reference data, specifically the moment tensor potential (MTP), combined with density functional theory (DFT) calculations and active learning loops that enable rapid sampling of MD trajectories. For SEI relevant materials (e.g., Li₂CO₃, bulk Li, LiPF₆, and Li₂EDC), our trained MTP models accurately capture the key structural properties (e.g., lattice parameters, elastic constants, or phonon spectra). For the dynamical properties of monoclinic Li₂CO₃ and amorphous Li₂EDC, the models are validated against previous theoretical predictions in the literature. Particularly, we illustrate the finite temperature effects on computing energy barriers. The determined mechanism of dominant diffusion carriers (Li vacancy, interstitial Li, and Li Frenkel pair) in Li₂CO₃ is highly consistent with DFT calculations. Furthermore, we show that the generated training datasets can be applied to train graph-neural-network (GNN)-based interatomic potentials that can further improve accuracy. The developed machine learning workflow provides a scalable approach for SEI modelling, enabling simulations at larger time and length scales to tackle the limitations of conventional DFT methods.

---

New concepts

The present study introduces a machine learning-driven framework for modelling the complex structure and dynamics of the solid electrolyte interphase (SEI) in lithium-ion batteries. By leveraging moment tensor potentials (MTPs) trained on density functional theory (DFT) data, the approach enables accurate and scalable molecular dynamics simulations. A key innovation lies in the integration of active learning loops and amorphous configurations, which iteratively refine the training set to capture a broader range of atomic environments efficiently. The workflow also demonstrates the transferability of generated datasets to graph neural network (GNN)-based force fields, offering a pathway to further improve predictive accuracy for SEI materials.

1. Introduction

Over the last decade, the widespread adoption of electric vehicles and portable electronics has surged dramatically, driven by the advancements in rechargeable Li-ion batteries (LIBs). However, further improvements are needed not only to enhance power capacity but also to improve safety, reliability, and lifespan across diverse applications.^1,2 One of the key issues in LIBs is the formation of kinetically stable solid electrolyte interphases (SEIs), which play a crucial role in battery degradation that results in capacity loss and increased resistance. Despite extensive research, the structure and evolution of the SEI remain elusive due to its complexity and the wide range of physical and chemical processes involved.³ The SEI forms through the sacrificial decomposition of the electrolyte and additive molecules during the initial charging of LIBs.⁴ The structure of the SEI is generally considered to consist of two distinct layers: an inner layer (close to the anode) primarily composed of fully reduced crystalline compounds (e.g., Li₂CO₃, Li₂O, and LiF) and an outer layer (close to the liquid electrolyte) composed of porous and non-crystalline phases such as organic dilithium ethylene dicarbonate (Li₂EDC) and dilithium butylene dicarbonate (Li₂BDC). In practice, the composition of the SEI is even more complex and influenced by various factors such as the choice of electrode materials, electrolyte compositions, salt concentrations, and reaction conditions such as the temperature and applied electric field.¹ A wide range of experimental techniques have been employed to characterize the SEI in detail.^5,6 However, gaining fundamental insights into SEI formation and growth remains challenging due to the inherent limitations of studying multiscale phenomena experimentally.

In battery materials research, molecular dynamics (MD) modelling has become a crucial tool for developing next-generation battery technologies.^7,8 Ab initio molecular dynamics (AIMD) based on density functional theory (DFT) has been widely used to model battery materials at the atomic scale, typically involving a few hundred atoms and picosecond time scales.⁹ However, AIMD struggles to model the SEI realistically due to its high computational cost and limited simulation scale.¹⁰ Empirically fitted force fields provide a complementary approach that can enable large-scale MD simulations by using analytic formulas to describe the interatomic interactions.^11–14 These models significantly reduce the computational cost as compared with DFT. However, the development of traditional force fields, which requires reparameterization of data points (experimental and/or quantum chemical data) with significant human involvement, is a tedious process for complex chemistries (particularly those involving bond breaking and formation), leading to a low transferability.

Machine learning-based force fields (MLFFs) are transforming the way we model complex material systems, offering an efficient yet highly accurate alternative to traditional interatomic potentials.^15–17 Recent progress in MLFF models for solid electrolyte materials in lithium-ion batteries has demonstrated their utility in probing structural stability, dopant effects, diffusivity, and activation energies across a wide range of chemistries.^18,19 By learning energy landscapes directly from first-principles calculations, MLFFs can bridge the gap between computational efficiency and quantum-level precision, making them particularly well-suited for studying intricate environments. Unlike classical force fields, which rely on predefined functional forms, MLFFs use flexible, data-driven models to capture a broad spectrum of atomic interactions. This adaptability is crucial for SEI materials, where transport properties depend on subtle variations in the chemical composition, structural disorder, and interfacial effects. Despite their promise, applying MLFFs to SEI systems remains an open challenge. One of the primary hurdles is the need for extensive, high-quality training datasets that comprehensively represent the diverse atomic configurations found in SEI materials. The SEI consists of a combination of organic and inorganic phases, often exhibiting both crystalline and amorphous regions, as well as dynamically evolving interfacial structures. A robust MLFF must accurately describe interactions among key elements such as Li, O, C, H, and F, while also capturing the complex behavior of molecular species like CO₃²⁻ and EDC²⁻ anions.

A major challenge in MLFF development is the computational cost associated with generating training data. First-principles calculations, particularly for large systems, are resource-intensive, limiting the feasible size of individual training structures to a few hundred atoms. Consequently, MLFF training datasets must be carefully curated to maximize coverage of relevant atomic interactions within these constraints. In practical applications, the architecture and technical setup of MLFFs—including the choice of reference data, the method of sampling, and the level of theory used—play a crucial role in determining the reliability and range of applicability of the final model. Several strategies have been proposed to assemble reference datasets from DFT or ab initio calculations, each with distinct advantages and limitations. For example, (i) AIMD sampling: running long molecular dynamics trajectories at finite temperatures to systematically explore relevant atomic configurations;^20,21 (ii) multiscale training: using lower-accuracy methods to generate an extensive dataset, followed by selective high-accuracy DFT refinements to enhance model precision;^20,22 (iii) active learning: dynamically improving the force field by incorporating novel configurations encountered during simulations;^23–26 (iv) enhanced sampling techniques: using biasing methods to accelerate the exploration of rare but important configurational states;²¹ and (v) perturbation-based methods: generating a diverse training set by displacing atomistic structures to capture high-energy states relevant for diffusion and reaction pathways.^17,27 Although these techniques have been successfully applied to various material systems, their optimal combination for SEI modelling remains unclear. The lack of standardized methodologies for constructing SEI-specific training datasets has hindered the development of MLFFs capable of capturing the full range of SEI phenomena. Additionally, validating MLFF predictions against DFT/experimental data remains an essential but underexplored aspect of this field. As computational resources continue to advance, MLFFs are poised to play a pivotal role in accelerating SEI research. However, achieving their full potential requires further innovation in dataset construction, model transferability, and validation frameworks. Addressing these challenges will be essential for leveraging MLFFs to gain deeper insights into the fundamental processes governing SEI formation, evolution, and ionic transport.

In this work, we demonstrate the utility of machine learning interatomic potentials (MLIPs) developed through systematically sampling amorphous structures via a batch of short MD simulations, coupled with an active learning scheme and DFT calculations, to efficiently and reliably model the major components of the SEI in LIBs. Direct comparison and validation of MLFF predictions against DFT benchmarks for key properties such as Li⁺ transport barriers, lattice stability, phonon spectra, and elastic moduli establish trust in the MLFF's predictive power for SEI chemistry. While the study builds upon existing machine learning frameworks such as MTP and MACE, its novelty lies in developing a domain-specific MLFF workflow tailored for the SEI—a chemically complex, structurally diverse, and dynamically evolving region that is poorly described by classical force fields.

2. Methodology

2.1 Setting up starting configurations

Traditionally, AIMD simulations constitute a straightforward way to explore configurational space starting with a selected structure. The temperature of simulations determines which regions of the potential energy surface (PES) and what energy ranges (according to the Boltzmann distribution) are explored in a long-time trajectory. However, for Li₂EDC and Li₂CO₃, exploration of the region relevant to rotation and translation of fluxional groups (EDC²⁻ and CO₃²⁻) cannot be achieved simply due to sluggish EDC²⁻ and CO₃²⁻ species, which usually need a very long-time simulation at elevated temperatures. On the other hand, the C–H, C–C and C–O bonds in EDC²⁻ might easily break at a temperature >120 °C according to previous experimental studies, leading to failure of sampling EDC²⁻ configurations.²⁸ Because of both points, different starting amorphous configurations with short MD simulations are used to explore various intramolecular interactions (e.g., EDC²⁻) and fluxional groups (CO₃²⁻). In addition, the machine learning force field parameters are fit for the configurational space of systems, bulk Li and LiPF₆.

The SEI exhibits a highly complex chemical environment, encompassing diverse compositions, interfaces, and stoichiometries. Consequently, its structural properties cannot be adequately represented by well-defined crystalline phases alone. In many cases, the interfacial regions of the SEI are inherently disordered at the atomic scale. To address this, we adopted a strategy that captures interactions between different compound types and overcomes the limitations of training machine learning potentials solely on near-crystalline configurations. Specifically, the initial structures were generated by randomly blending and packing varying amounts of Li₂EDC, Li₂CO₃, LiPF₆ and Li species into a supercell using an initialisation scheme proposed by Vilhelmsen and Hammer, Fig. S1.²⁹ The supercell involves atoms with a maximum number of 8 Li₂EDC, 16 Li₂CO₃, 64 Li molecules, etc. Initially, each system was set up at low density with a large supercell (Fig. S1). Then it was condensed over 5 ps at 10 K with a high pressure of 50 kbar in the NPT ensemble, and an optimization was conducted to relax the atomic positions and internal stress. These calculations were performed using a universal machine learning force field implemented in MACE.³⁰ In this work, 120 starting structures were generated this way to evaluate the training set via the loops described in the following section (Scheme 1). The results and discussion section demonstrates that the potential exhibits strong predictive capability for various properties near equilibrium states (Fig. 1).

Scheme 1 Active learning bootstrapping iterations (i)–(iv) with combinations of density functional theory (DFT), molecular dynamics (MD) and moment tensor potential (MTP). The critical extrapolation grades (γ) are set to 10 for terminating MD simulations and 2 for selecting configurations to extend the training set.

Fig. 1 Typical structure of the SEI with inorganic and organic components, e.g., Li₂CO₃ and Li₂EDC. Brown – carbon, red – oxygen, white – hydrogen, silver – fluorine, and light purple – phosphorus.

2.2 MLFF training procedure

A strategy for sufficiently sampling desired configurations of atomistic structures is outlined in Scheme 1. The moment tensor potential (MTP) framework²⁶ and a bootstrapping technique implemented in the MLIP-3 package³¹ are used to assemble the training set via sampling configurations directly from a targeted atomistic simulation. This concurrent learning procedure allows us to terminate the simulation according to extrapolation criteria,³¹ retrain a potential, and restart a simulation. Given an initial structure, the active learning loop may require several iterations to produce a robust machine learning potential. Each iteration of the loop typically follows five key steps, as illustrated in Scheme 1: (i) initialization: a preliminary training dataset is generated from a short AIMD simulation, which is used to construct a pre-trained potential. (ii) Exploratory simulation: molecular dynamics is performed using the current potential, during which extrapolative configurations—those that fall outside the model's reliable prediction space—are actively identified. The MTP-3 software^26,31 employs the extrapolation grade γ, derived from the D-optimality criterion,^32–34 as a quantitative measure of how ‘unfamiliar’ a new atomic configuration is relative to those already in the training set. If the MD simulation terminates upon exceeding the predefined extrapolation threshold, γ > 10 commonly used for MTP training,^26,31 an update to the active training set is required. To accomplish this, the Maxvol algorithm is employed to select the most informative configurations from the terminated trajectory.³¹ In the Li₂CO₃ system, test runs with γ = 8 and γ = 10 yielded nearly identical results, with differences of only 9 meV Å⁻¹ in force RMSE, 0.2 meV per atom in energy RMSE, and 10 structures in training set sizes. (iii) DFT calculations: the configurations with γ > 2 are selected and evaluated using the DFT single point calculation to obtain accurate reference data including energies, forces, and stresses. These configurations with DFT information are then added to the training dataset. (iv) Model update: the MTP is retrained with the expanded dataset, and the next iteration using newly trained MTP begins from step (ii). Each cycle incrementally broadens the model's training domain and enhances the stability and transferability of the potential—allowing longer and more reliable simulations without triggering extrapolation failures. This iterative loop continues until the simulation completes without surpassing the critical extrapolation grade γ_break.²⁶ Notably, for a fixed thermodynamic condition, each iteration may initiate a different trajectory, as the molecular dynamics are typically restarted with a new random velocity seed.

In the MTP framework, the potential energy of an atomic structure is approximated as the sum of site energies for individual atoms, where each site's energy is expanded linearly using a set of basis functions. MTP introduces degree-like measures, referred to as “levels,” to determine which basis functions are included in the interatomic potentials. In our implementation, we employ MTP at “level 16,” corresponding to a mid-sized basis set, which enables large-scale MD simulations with sufficient efficiency. This is particularly advantageous for investigating SEI materials, where computational speed is critical. We have observed that this strategy can result in sufficient sampling for intended applications of the final MTP model involving MD simulations for equilibrium or close to equilibrium properties (e.g., radial distribution functions and lithium-ion diffusion coefficients).

Starting with 120 initial configurations generated via the scheme discussed in Section 2.1, the active learning process (see Section 2.2) ultimately produced the final training set, comprising 8686 structures, which enabled the generation of 50-ps NVT and NPT trajectories at various temperatures below 1000 K and pressures of 1 bar and 100 kbar. Training at MTP “level 16” yields a root-mean-squared error (RMSE) of 8 meV per atom for energy and a RMSE of 0.34 eV Å⁻¹ for forces.

2.3 DFT and molecular dynamics methods

The DFT calculations were carried out using the Vienna ab initio Simulation (VASP) package.^35,36 The projected augmented wave method and the generalized gradient approximation of Perdew–Burke–Ernzerhof (PBE) type for the exchange correlation functional are employed.³⁷ An energy cutoff of 520 eV and an automatic k-mesh generation are used to converge the energy and atomic force. The number of k-points in the direction of the first, second and third reciprocal lattice vectors is determined by setting the smallest allowed spacing (0.35 Å⁻¹) between k-points.

All of the MD simulations were performed using the MTP force field³¹ implemented as a calculator in the LAMMPS package³⁸ using a constant time step of 0.5 fs. All simulations were performed using NPT or NVT ensembles, as appropriate, using a Nosé/Hoover temperature thermostat and a Nosé/Hoover pressure barostat, including an ensemble of velocities from a random number generator with the specified seed at the specified temperature. The relaxation time constants were chosen to be 100 fs for the thermostat and 500 fs for the barostat. In the production stage of the MD simulations, the initial geometries of the amorphous Li₂EDC were created using the disordered system builder as implemented in ASE proposed by Vilhelmsen and Hammer²⁹ (see Fig. S1 and S6 for more details).

2.4 MLIP performance on test sets and validation using DFT

Starting with 120 initial configurations (see Section 2.1), we perform 50-ps NPT MD simulations at 300 K and 1 bar for each structure using the trained MTP model. We extract a test set of 20 000 sampled structures from those trajectories to evaluate the model performance. We assess the model's accuracy by comparing the absolute energies (Fig. 2a) and force components (Fig. 2b) against the reference level from theory (DFT with the PBE functional). The test root-mean-squared errors (RMSEs) for energy per atom and force components are 9 meV and 0.21 eV Å⁻¹, respectively—closely matching the RMSEs of the training set. This indicates an accurate and robust model is achieved. More importantly, nevertheless, the more stringent and useful tests for our model are quantitatively accurate predictions of various properties (e.g., lattice parameters, elastic constants, phonon spectra, radial distribution functions and lithium-ion diffusion coefficients) of materials compared to DFT results and previous theoretical predictions in the literature. These tests are discussed in Section 3.

Fig. 2 Correlation of test set energies (a) and atomic force x components (b) of the DFT level with the MTP predictions.

3. Results and discussion

3.1 MLFF performance in predicting material properties at local energy minima

We present validation calculations for the lattice parameters and elastic constants of various inorganic and organic materials, including monoclinic Li₂CO₃, bcc bulk Li, hexagonal LiPF₆, and organic Li₂EDC, in Fig. 3. Additionally, we consider two high-pressure phases: hexagonal Li₂CO₃ and fcc bulk Li. Experimental data for the atomic structure of crystalline Li₂EDC are unavailable; however, XRD measurements suggest an orthorhombic structure with lattice constants (a = 14.10 Å, b = 12.66 Å, c = 5.20 Å), accommodating six Li₂EDC molecules based on the estimated density.²⁸ Since our goal is to compare MLFF predictions with DFT-derived lattice parameters, we limit our calculations to a unit cell containing a single Li₂EDC molecule (Fig. 3). The MTP calculations accurately reproduce the DFT-derived lattice parameters, with absolute deviations within 2% (Table S1). Moreover, elastic constants are fundamental material properties that are closely linked to the second derivative of the potential energy and chemical bonds between atoms. The MLFF model provides reasonable predictions of the elastic properties across different materials, closely matching the DFT values and effectively capturing variations among them (Table S2). Overall, the MLFF model demonstrates good agreement with DFT calculations of material properties at local energy minima. Notably, the training database does not explicitly contain configurations from these ordered phases, underscoring the model's accuracy and transferability.

Fig. 3 Atomic structure of monoclinic Li₂CO₃ (a), hexagonal Li₂CO₃ (b), organic Li₂EDC (c), bcc bulk Li (d), fcc bulk Li (e), and hexagonal LiPF₆ (f). Red – oxygen, green – lithium, brown – carbon, white – hydrogen, light purple – phosphorus, and silver – fluorine. Solid lines indicate the unit cells.

3.2 MLFF performance in predicting dynamic properties

Simulations on dynamic properties for the SEI require large atomic models typically containing thousands of atoms. Therefore, we tested whether the machine-learned potential, trained on small cells, can scale to larger models and extended time scales while retaining near-DFT accuracy. To verify this, we benchmarked the potential on a series of progressively larger Li₂EDC and Li₂CO₃ supercells. The models are constructed following the scheme in Fig. S1. As an example, we tracked a long MD trajectory of a large supercell containing amorphous Li₂EDC (Fig. 4a). Throughout the MD trajectory, the MTP energies track the DFT values within roughly 10 meV per atom. As shown in Fig. 4b and c, the energy RMSE is essentially unchanged across cell sizes, while the force RMSE increases only slightly, from 0.34 to 0.43 eV Å⁻¹, for larger supercells. These results indicate that the MTP remains stable and gives accurate predictions over long MD simulations on large and previously unseen configurations.

Fig. 4 Correlation of test set energies and force at the DFT level with the MTP predictions. (a) Atomic structure of an amorphous Li₂EDC in a large cell (left), and the total energy per atom calculated using DFT and MTP from a long MD trajectory. Energy (black) and force (blue) root-mean-squared errors in (b) Li₂EDC and (c) Li₂CO₃ for a series of supercells with different atomic sizes, with the horizontal axes indicating the different number of atoms in supercells.

3.2.1 Li diffusivity within Li₂EDC and Li₂CO₃. Previous theoretical studies using DFT calculations have proposed three primary mechanisms for Li transport in SEI components:^9,39,40 (i) Li hopping between neighboring sites,⁴¹ (ii) Li hopping assisted by Li vacancies,⁴⁰ and (iii) interstitial Li transport via either a “knock-off” mechanism⁹ or open-channel diffusion.³⁹ The predominant Li migration mechanism in the simulations of stoichiometric Li₂EDC and Li₂CO₃ primarily involves Li hopping between neighboring sites. Our MD simulations show that, a Li atom adjacent to a vacancy (created when another Li moves to a different position, Fig. S9) tends to fill the vacant position during Li migration. Due to the strong interaction between Li⁺ and CO₃²⁻, stoichiometric Li₂CO₃ exhibits a high activation energy hop barrier of approximately 1.2 eV (Fig. S8). Nudged elastic band (NEB) calculations using MTP further indicate that the formation of a neutral Li Frenkel pair via Li migration into an open channel along the [010] direction requires an energy barrier of 1.19 eV (Fig. 5a). Unlike pure Li₂CO₃, the presence of a Li vacancy and interstitial Li significantly enhances diffusion, particularly at lower temperatures, through the hopping mechanism assisted by the Li vacancy and the knock-off diffusion mechanism, respectively. As shown in Fig. 5b and c, the energy barriers calculated using MLFFs and DFT via the NEB method closely align with the values estimated from diffusivity. The barriers also align with the previous DFT-determined range of 0.20 to 0.60 eV for interstitial Li transport,^9,39,42 demonstrating the reliability of our trained MLFFs in capturing SEI properties.

Fig. 5 Energy profile of the pathway following (a) the generation of a neutral Li Frenkel pair, (b) the neutral Li vacancy diffusion, and (c) the neutral Li diffusion via a knock-off mechanism. Black – DFT, blue – MTP, and yellow – MACE. The nudged elastic band method⁴³ as implemented in ASE⁴⁴ was used to locate the energy barrier of Li diffusion in a 1 × 2 × 2 supercell of Li₂CO₃. The atomistic structures of initial, transition, and final states are shown in Fig. S10.

Finite-temperature dynamics offer critical insights into lithium-ion transport within the SEI materials, which is an important method for understanding batteries based on electrode–electrolyte systems. The experimentally determined activation energy (∼0.60 eV) near the anode is expected to be an average value over the multiphase SEI including different interfaces.^42,45 Accurately predicting lithium-ion transport across different phases requires well-trained force fields. As an example, we consider models consisting of 1440 atoms from a 3 × 5 × 4 supercell of monoclinic Li₂CO₃, as well as an amorphous model containing 72 Li₂EDC molecules (Fig. S6 in the SI) with 1152 atoms. MD simulations were performed at elevated temperatures (600–900 K for Li₂CO₃ and 400–600 K for Li₂EDC) to ensure that the system reached a diffusive regime, allowing the analysis of temperature-dependent transport behaviours. The self-diffusion coefficient (D) was calculated using the Einstein relation at temperatures where diffusion was observed in NPT simulations at a pressure of 1 bar:

Here, MSD(t) represents the mean-square displacement of Li over time t, with 〈·〉 denoting the ensemble average. The parameter n corresponds to the dimensionality of the system, where n = 3 for three-dimensional diffusion in our model. Note that, at low temperatures (T < 400 K for Li₂EDC and T < 600 K for Li₂CO₃), atomic vibrations dominate, with few diffusion hopping events observed within the simulation timescale. Before evaluating the mean-square displacement of Li, the atomic models were equilibrated for at least 80 ns for the amorphous Li₂EDC supercell (Fig. S6) and 1 ns for the crystalline Li₂CO₃ supercell.

The diffusive behaviour of Li ions is evident from the linear increase of MSD with time (Fig. S7), and the temperature dependence follows an Arrhenius-like relationship, appearing linear in the plot of log(D) vs. 1/T (Fig. 6). While the MSD is higher at elevated temperatures, Li hopping events are significantly reduced at room temperature, leading to much lower diffusivity. At 300 K, our simulations yield lithium-ion diffusivities of 1.25 × 10⁻¹⁵ m² s⁻¹ for Li₂EDC and 1.32 × 10⁻²³ m² s⁻¹ for Li₂CO₃ (Fig. 6). The diffusion coefficient in Li₂EDC is consistent with the value (∼10⁻¹⁶ m² s⁻¹) obtained using the revised many-body polarizable APPLE&P force field.¹² To the best of our knowledge, only Tasaki et al.⁴⁶ have previously reported Li diffusion in crystalline Li₂CO₃ at room temperature, using the COMPASS force field with parameters by defined Accelrys, Inc. In detail, they used a Li₂CO₃ crystal model with a simulation cell of approximately 20 Å per side and performed 1 ns NPT ensemble simulations. In our study, we employed a larger supercell (∼25 Å) and extended the simulation time to 20 ns under the NPT ensemble (Fig. S7) in order to evaluate the Li⁺ diffusion coefficient in Li₂CO₃. However, our calculated Li diffusion coefficient for Li₂CO₃ is approximately eight orders of magnitude lower than the value, 9 × 10⁻¹⁵ m² s⁻¹, reported by Tasaki et al.⁴⁶ This discrepancy may arise from differences in computational models or the methodologies used to estimate diffusivity. Our calculated value (1.32 × 10⁻²³ m² s⁻¹ for Li₂CO₃) was used to derive the activation energy barrier (Fig. S8a), which was further validated via DFT calculations presented in Fig. 5. Moreover, a comprehensive study of Li transport in Li₂EDC, combining experimental measurements with the APPLE&P force field,¹² has indicated that the diffusion coefficient of 8 × 10⁻¹² m² s⁻¹ reported by Tasaki et al. is significantly overestimated compared to our more consistent value of ∼10⁻¹⁵ m² s⁻¹.

Fig. 6 Diffusion coefficients of Li obtained from MD simulations of pure Li₂CO₃ (a), Li₂CO₃ with interstitial Li (b), Li₂CO₃ with the Li vacancy (c) and amorphous Li₂EDC (d). The dashed lines highlight the implied Arrhenius-like dependence for Li transport.

The obtained activation energy barrier for Li-ion diffusion is 0.64 eV for Li₂EDC (Fig. S8), which closely aligns with the previously reported value of 0.66 eV using the classical force field.¹² For crystalline Li₂CO₃, the barrier is 1.19 eV, which is highly consistent with our DFT prediction and falls within the range of 0.78–1.34 eV measured using solid-state NMR spectroscopy.⁴⁷ Within the voltage range of lithium-ion batteries, excess interstitial Li is expected to dominate the diffusion process due to its lower formation energy compared to other defects (e.g., Li vacancy and Li Frenkel pair),⁹ making it a key factor in estimating energy barriers.^9,39,40 Our simulations reveal a more rapid linear increase in the MSD curve (Fig. S7 in the SI) in the presence of an extra interstitial Li in the Li₂CO₃ model with 1440 atoms, confirming its impact on diffusivity. The concentration of interstitial Li-doped Li₂CO₃ is about 7.1 × 10²⁵ m⁻³, which is close to the concentration used in SEI growth modelling.⁴⁸ Fig. 6b further illustrates the Arrhenius-like relationship. The diffusivity of the interstitial Li model at 300 K is 4.25 × 10⁻¹³ m² s⁻¹, comparable to the values obtained using AIMD (7.6 × 10⁻¹⁵ m² s⁻¹)⁴² and the standard 6–12 Lennard-Jones (LJ) potential (3.3 × 10⁻¹⁶ m² s⁻¹).⁴¹ The activation energy for interstitial Li diffusion in Li₂CO₃ is found to be 0.18 eV. Additionally, Li vacancies are frequently reported as point defects that facilitate Li transport in Li₂CO₃.⁴⁹ At 300 K, the Li vacancy diffusivity is 8.0 × 10⁻¹⁴ m² s⁻¹ with an activation energy of 0.24 eV, with both values comparable to those of interstitial Li. The close agreement between energy barriers estimated from DFT-based NEB calculations and those obtained from MLFF-based MD simulations further validates the high accuracy of our trained model.

The diffusivities of interstitial Li and Li vacancies in Li₂EDC at 300 K are 1.65 × 10⁻¹⁴ m² s⁻¹ and 6.11 × 10⁻¹⁵ m² s⁻¹, respectively, with the corresponding energy barriers of 0.49 eV and 0.56 eV. These values suggest that the organic Li₂EDC exhibits lower ionic mobility than the inorganic Li₂CO₃, aligning with experimental findings where the stripping of the Li⁺ solvation sheath⁵⁰ and the diffusion of Li⁺ from the organic to the inorganic layer⁹ are identified as key steps in anode processes.

We explicitly calculated the ionic conductivities of both Li₂CO₃ and Li₂EDC using the Nernst–Einstein relation: σ = nz²q²D/k_BT, where σ is the ionic conductivity, n is the charge carrier concentration, D is the diffusion coefficient, T is the temperature, z = 1 is the charge number for Li⁺, q = 1.602 × 10⁻¹⁹ C is the elementary charge, and k_B = 1.381 × 10⁻²³ J K⁻¹ is the Boltzmann constant. We adopt a larger supercell containing more atoms for the ionic conductivity calculations so that the system size is more comparable to other theoretical work. Previous studies^51–56 have shown that ionic conduction in inorganic Li₂CO₃ primarily occurs via interstitial Li defects.

The charge carrier concentration n of Li interstitials in Li₂CO₃ computed using a supercell including 3457 atoms is 2.98 × 10²⁵ m⁻³. The calculated ionic conductivity based on the computed diffusion coefficients (1.76 × 10⁻¹³ m² s⁻¹ at 300 K) is 3.25 × 10⁻⁵ S m⁻¹. Available experimental studies reported that the ion conductivity of Li₂CO₃ is ∼10⁻⁸ S m⁻¹.^51–53 Nevertheless, the theoretical reported values based on DFT and the continuum approach typically range from 10⁻⁶ to 10⁻⁸ S m⁻¹.^54–56 For Li₂EDC, unlike Li₂CO₃, the dominated ion conduction mechanism remains experimentally unclear. Therefore, both intrinsic Li and Li interstitials might contribute to the conduction. With this assumption, the estimated charge carrier concentrations for Li₂EDC using a supercell including 2688 atoms range from 3.50 × 10²⁵ to 1.18 × 10²⁸ m⁻³. The calculated ionic conductivities range from 7.60 × 10⁻⁴ to 2.26 × 10⁻⁶ S m⁻¹ based on the computed diffusion coefficients (1.04 × 10⁻¹⁴ m² s⁻¹ at 300 K). The experimental ionic conductivity is ∼10⁻⁷ S m⁻¹, while a forcefield based modelling gives values ranging from 10⁻⁶ to 10⁻⁸ S m⁻¹.¹²

Despite the calculated value of Li₂CO₃ being higher than the experimental value, our result is in good agreement with other theoretical work, whereas the Li₂EDC result agrees with both experimental and other theoretical works. We note that the ionic conductivity calculated from modelling in the literature tends to be higher. This discrepancy could primarily stem from the fact that it is computationally challenging to model the experimental charge carrier's concentration which requires a very large supercell. Therefore, a higher concentration of charge carrier is assumed in modelling (i.e., adopt a relatively smaller supercell), which could result in higher ionic conductivity. Nevertheless, the trend obtained from the present work remains consistent with other modelling and experimental work.

3.2.2 Li transport between Li(111) and Li₂CO₃(001) interfaces. We present a case study of Li(110)/Li₂CO₃(001) and Li(111)/Li₂CO₃(001) interfaces (Fig. S12) to demonstrate the applicability and extensibility of our MLFF workflow to interfacial systems. The interface energies are 0.44 and 0.35 J m⁻² calculated using the MTP model for both interfaces, respectively, which are consistent with DFT values, 0.52 J m⁻² and 0.46 J m⁻², respectively. We have observed that the Li(111)/Li₂CO₃(001) interface is denser than the Li(110)/Li₂CO₃(001) interface due to Li atoms from metal Li form more bonds with the O of CO₃ at Li₂CO₃(001), leading to a more stable interface between Li(111) and Li₂CO₃(001) surfaces.

In addition, previous DFT studies suggest that the Li(001)/Li₂CO₃(001) interface exhibits relatively high mechanical strength with a reported interface energy of 0.49 J m⁻²,⁵⁸ yet Li⁺ migration across it has not been thoroughly investigated. Our MD simulation for Li(001)/Li₂CO₃(001) interface shows that Li(001) changes to Li(111) after 5 ns at 300 K, indicating that Li(111) is the lowest energy surface. Therefore, we have investigated Li-ion transport across the Li(111)/Li₂CO₃(001) interface (Fig. 7). MLFF-MD simulations performed under an applied electric field show that Li transport takes place from Li(111) to Li₂CO₃(001), as shown in Fig. 7. With the application of the electric field (right to left), the migration process takes place along the atomic chain (the five Li atoms labelled by numbers 1 to 5 in Fig. 7). Like the diffusion pathway in the “knock-off mechanism”, these Li atoms at sites from 5 to 1 take turns as the diffusing interstitials, move in four knock-off steps, and finally push the Li atom to the left side of Li₂CO₃.

Fig. 7 MLFF-MD simulations performed at 300 K under an applied electric field (0.03 eV e⁻¹ Å⁻¹) using the QEq method⁵⁷ to reveal lithium transport from Li(111) to Li₂CO₃(001). (a) Atomic structure of the Li(111)/Li₂CO₃(001) interface. Snapshots of the MD trajectory at (b) 1.4 ps and (c) 1.5 ps illustrate the migration process.

3.2.3 Li-related radial distribution functions. We further analyze the structures of amorphous Li₂EDC and crystalline Li₂CO₃ by examining their radial distribution functions (RDFs) obtained from MD simulations at 300 K and elevated temperatures of 500 K (Li₂EDC) and 700 K (Li₂CO₃). As shown in Fig. 8, the first Li–Oc peak for both materials appears at around 1.96 Å, indicating that the Li–Oc distance is primarily governed by electrostatic interactions between the Li cation and the carbonyl oxygen atom. Each Li ion coordinates with approximately four oxygen atoms, bridging aggregates of EDC²⁻ and CO₃²⁻ anions, as illustrated in Fig. 8b and d. The first Li–Li RDF peak shows that Li atoms are more widely spaced in Li₂CO₃ (3.00 Å, Fig. 8c) compared to Li₂EDC (2.89 Å, Fig. 8a). At elevated temperatures, the probability of Li coordinating with oxygen decreases slightly as compared to that at 300 K, reflecting increased Li mobility. In Li₂EDC, Li atoms exhibit a strong preference for coordinating with carbonyl oxygen (Oc) rather than ether oxygen (Oe), as confirmed by the Li coordination number in Fig. 8b. The Li–Oe peak appears at a slightly larger distance (2.12 Å) compared to the Li–Oc peak. The extended Li coordination shell, defined by the minimum following the first Li–Li RDF peak at 3.91 Å, indicates that Li⁺ remains surrounded by Oc in this secondary coordination environment. These findings are in strong agreement with previous simulations.^11,12 Additionally, due to the low interstitial-to-lattice Li ratio (1/480), the RDFs of pure and interstitial-doped Li₂CO₃ remain indistinguishable (Fig. S11).

Fig. 8 Radial distribution functions (a) and (c) of Li–X (X = Oc, Oe and Li) and the coordination numbers of Li for Li–X (b) and (d) obtained from MD simulations of amorphous Li₂EDC and crystalline Li₂CO₃ at different temperatures.

3.2.4 Improvement of the MTP model. Training at MTP “level 16” (92 descriptors) yields excellent performance for energy predictions, with a root-mean-squared error (RMSE) of 8 meV per atom. However, the RMSE of force prediction is 0.34 eV Å⁻¹. On a 128-core system (Dual-CPU AMD EPYC 7713), this setup achieves a physical simulation time of approximately 10 nanoseconds per day for a supercell containing ∼1000 atoms. Higher-order tensor terms only slightly improve accuracy but at a super-linear increase in computational cost. For instance, increasing the MTP level to 22 (288 descriptors) reduces the RMSE to 6 meV per atom for energy and 0.31 eV Å⁻¹ for forces, but at the cost of tripling the computational time. Given the chemically complex environments within SEI materials, MTP—which relies on a fixed basis expansion—performs well for smooth, interpolative energy landscapes but may struggle with capturing non-polynomial energy variations. As shown in Fig. 5, MTP exhibits energy barrier deviations of 0.03 to 0.07 eV for Li diffusion in Li₂CO₃ when compared to DFT calculations.

We have therefore further developed and validated an equivariant neural network model named MACE,³⁰ to assess whether the training dataset generated through iterative loops based on MTP models can be leveraged to capture subtle interatomic interactions with higher fidelity. The goal is to improve transferability across diverse material environments while enhancing accuracy. Training at MACE using about one million parameters yields a RMSE of 2 meV per atom for energy and a RMSE of 0.07 eV Å⁻¹ for forces (Fig. S2). Compared to the MTP model, the MACE model also achieves significantly lower RMSEs for both energy and forces on test sets (Fig. S2 and S3). Like the MTP model, the MACE model demonstrates excellent agreement with DFT calculations for lattice parameters (Table S1), elastic constants (Table S2), and phonon spectra (Fig. S3). Additionally, the MACE model improves predictions of energy barriers for Li transport via three mechanisms in Li₂CO₃ (Fig. 5) compared to the MTP model. However, despite these improvements, the computational cost of MACE remains a limiting factor. In large supercells (∼1440 atoms), our current implementation achieves only ∼180 picoseconds of MD simulation per day using one GPU CAR on a machine (Nvidia A100 40 GB). In contrast, MTP models indicate that equilibrating the amorphous Li₂EDC structure requires at least 60 nanoseconds of MD simulations (Fig. S6), highlighting a major challenge in applying MACE to SEI systems. This computational bottleneck restricts its feasibility for long-timescale simulations, such as studying ion transport across grain boundaries in SEI's organic or inorganic phases, which evolve due to dissolution, redeposition, and densification of different phases. Collectively, both methods highlight the importance of developing a diverse suite of MLFFs tailored to the specific accuracy-efficiency requirements of different battery components. Future extensions may benefit from hybrid strategies that combine the strengths of local-descriptor models like MTP with other more expressive approaches such as Sparse Gaussian Process Regression (SGPR). SGPR, as demonstrated in recent studies^59,60 on solid electrolyte materials, provides not only reliable accuracy but also intrinsic uncertainty quantification.

While our study focuses on Li₂CO₃ and Li₂EDC as representative SEI components due to their well-characterized roles in the inorganic and organic phases of the SEI, we acknowledge that the real SEI is a heterogeneous and dynamic mixture of various species, including LiF, Li₂O, organic polymers, and decomposition products including battery additives.^4–6,61 Our current framework is designed to be extensible, and the active learning workflow demonstrated here can be applied to incorporate additional SEI components in future studies. This work represents an important first step toward building more comprehensive machine learning-based models for the complex multicomponent nature of the SEI.

4. Conclusions

In this work, we have developed a machine learning interatomic potential approach to address the fundamental mechanism of Li transport in the SEI in lithium-ion batteries. We showed that the developed machine learning-based force fields are not only applicable to the major components of the SEI (i.e., Li₂CO₃ and Li₂EDC) but also generalizable to be able to include other components. e.g. LiF and Li₂O. The model trained based on active learning loops accurately reproduced the lattice parameters, elastic constants and phonon spectra across various SEI components, demonstrating its reliability and transferability. Through MD simulations, specifically focusing on amorphous Li₂EDC and crystalline Li₂CO₃, we revealed the lithium-ion diffusion mechanisms in different SEI phases. The obtained diffusivity values and activation energy barriers were in good agreement with experimental and computational benchmarks in the literature, further validating our approach. Overall, the work has demonstrated the potential of MLIPs in enabling large-scale high-accuracy modelling of the complex SEI structures in LIBs. By capturing the complex chemistry and ionic transport behaviour with near-DFT precision, our approach paves the way for predictive simulations that can guide the design of SEI structures for high-performance lithium-ion batteries.

Author contributions

W. Q. Li: conceptualization, analysis, visualization, and original draft. G. Wu, Y. H. Lau, and J. M. Arce-Ramos: analysis, review and editing. M. F. Ng: conceptualization, funding acquisition, review and editing, and supervision.

Conflicts of interest

The authors declare that they have no conflicts of interest.

Data availability

Relevant data (training data and MLIP models) are available in https://github.com/gitliwq/LiCOHPF_database_1.git.

Supplementary information: Tables and figures illustrating the structures and properties of various materials and systems. See DOI: https://doi.org/10.1039/d5mh01343g.

Acknowledgements

We acknowledge support from the A*STAR MTC programmatic funding (grant no. M23L9b0052) and the Singapore National Research Foundation's 2023 Singapore-China Joint Flagship Project (Clean Energy). We acknowledge the National Supercomputing Center (NSCC) Singapore and Agency for Science, Technology and Research (A*STAR) Computational Resource Centre (A*CRC) of Singapore for the use of its high-performance computing facilities. Wen-Qing thanks Dr Benjamin Chen for his technical support with the machine learning-based force fields.

References

1. J.-M. Tarascon and M. Armand, Nature, 2001, 414, 359–367.
2. J. B. Goodenough and K.-S. Park, J. Am. Chem. Soc., 2013, 135, 1167–1176.
3. Y. Peng, R. Tamate and K. Nishikawa, ChemElectroChem, 2024, 11, e202400278.
4. X. Xie, E. W. Clark Spotte-Smith, M. Wen, H. D. Patel, S. M. Blau and K. A. Persson, J. Am. Chem. Soc., 2021, 143, 13245–13258.
5. A. M. Tripathi, W.-N. Su and B. J. Hwang, Chem. Soc. Rev., 2018, 47, 736–851.
6. Y. Zhou, M. Su, X. Yu, Y. Zhang, J.-G. Wang, X. Ren, R. Cao, W. Xu, D. R. Baer, Y. Du, O. Borodin, Y. Wang, X.-L. Wang, K. Xu, Z. Xu, C. Wang and Z. Zhu, Nat. Nanotechnol., 2020, 15, 224–230.
7. A. Wang, S. Kadam, H. Li, S. Shi and Y. Qi, npj Comput. Mater., 2018, 4, 15.
8. D. Diddens, W. A. Appiah, Y. Mabrouk, A. Heuer, T. Vegge and A. Bhowmik, Adv. Mater. Interfaces, 2022, 9, 2101734.
9. S. Shi, P. Lu, Z. Liu, Y. Qi, L. G. Hector, H. Li and S. J. Harris, J. Am. Chem. Soc., 2012, 134, 15476–15487.
10. Y. Li, K. Leung and Y. Qi, Acc. Chem. Res., 2016, 49, 2363–2370.
11. D. Bedrov, O. Borodin and J. B. Hooper, J. Phys. Chem. C, 2017, 121, 16098–16109.
12. O. Borodin, G. V. Zhuang, P. N. Ross and K. Xu, J. Phys. Chem. C, 2013, 117, 16749.
13. N. Takenaka, Y. Suzuki, H. Sakai and M. Nagaoka, J. Phys. Chem. C, 2014, 118, 10874–10882.
14. L. Raguette and R. Jorn, J. Phys. Chem. C, 2018, 122, 3219–3232.
15. Y. Shao, L. Knijff, F. M. Dietrich, K. Hermansson and C. Zhang, Batteries Supercaps, 2021, 4, 585–595.
16. V. L. Deringer, JPhys: Energy, 2020, 2, 041003.
17. S. Dajnowicz, G. Agarwal, J. M. Stevenson, L. D. Jacobson, F. Ramezanghorbani, K. Leswing, R. A. Friesner, M. D. Halls and R. Abel, J. Phys. Chem. B, 2022, 126, 6271–6280.
18. M. Ha, A. Hajibabaei, D. Y. Kim, A. N. Singh, J. Yun, C. W. Myung and K. S. Kim, Adv. Energy Mater., 2022, 12, 1–9.
19. S. Y. Willow, A. Hajibabaei, M. Ha, D. C. Yang, C. W. Myung, S. K. Min, G. Lee and K. S. Kim, Chem. Phys. Rev., 2024, 5, 041307.
20. O. T. Unke, S. Chmiela, H. E. Sauceda, M. Gastegger, I. Poltavsky, K. T. Schütt, A. Tkatchenko and K. R. Müller, Chem. Rev., 2021, 121, 10142–10186.
21. D. Yoo, J. Jung, W. Jeong and S. Han, npj Comput. Mater., 2021, 7, 1–9.
22. G. Vazquez, D. Sauceda and R. Arróyave, Acta Mater., 2024, 276, 120137.
23. G. Csányi, T. Albaret, M. C. Payne and A. De Vita, Phys. Rev. Lett., 2004, 93, 175503.
24. R. Jinnouchi, J. Lahnsteiner, F. Karsai, G. Kresse and M. Bokdam, Phys. Rev. Lett., 2019, 122, 225701.
25. B. W. J. Chen, X. Zhang and J. Zhang, Chem. Sci., 2023, 14, 8338–8354.
26. I. S. Novikov, K. Gubaev, E. V. Podryabinkin and A. V. Shapeev, Mach. Learn. Sci. Technol., 2021, 2, 025002.
27. Y. Zhang, H. Wang, W. Chen, J. Zeng, L. Zhang, H. Wang and W. E, Comput. Phys. Commun., 2020, 253, 107206.
28. K. Xu, G. V. Zhuang, J. L. Allen, U. Lee, S. S. Zhang, P. N. Ross and T. R. Jow, J. Phys. Chem. B, 2006, 110, 7708–7719.
29. L. B. Vilhelmsen and B. Hammer, J. Chem. Phys., 2014, 141, 044711.
30. I. Batatia, D. P. Kovacs, G. N. C. Simm, C. Ortner and G. Csanyi, in Advances in Neural Information Processing Systems, ed. A. H. Oh, A. Agarwal, D. Belgrave and K. Cho, 2022.
31. E. Podryabinkin, K. Garifullin, A. Shapeev and I. Novikov, J. Chem. Phys., 2023, 159, 084112.
32. K. Gubaev, E. V. Podryabinkin and A. V. Shapeev, J. Chem. Phys., 2018, 148, 241727.
33. K. Gubaev, E. V. Podryabinkin, G. L. W. Hart and A. V. Shapeev, Comput. Mater. Sci., 2019, 156, 148–156.
34. E. V. Podryabinkin and A. V. Shapeev, Comput. Mater. Sci., 2017, 140, 171–180.
35. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 14251–14269.
36. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
37. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
38. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, Comput. Phys. Commun., 2022, 271, 108171.
39. H. Iddir and L. A. Curtiss, J. Phys. Chem. C, 2010, 114, 20903–20906.
40. Y. Chen, M. Huo, L. Song and Z. Sun, RSC Adv., 2014, 4, 42462–42466.
41. L. Benitez and J. M. Seminario, J. Electrochem. Soc., 2017, 164, E3159–E3170.
42. Z. Ahmad, V. Venturi, H. Hafiz and V. Viswanathan, J. Phys. Chem. C, 2021, 125, 11301–11309.
43. G. Henkelman and H. Jónsson, J. Chem. Phys., 2000, 113, 9978–9985.
44. A. Hjorth Larsen, J. Jørgen Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. Bjerre Jensen, J. Kermode, J. R. Kitchin, E. Leonhard Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. Bergmann Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiøtz, O. Schütt, M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. Zeng and K. W. Jacobsen, J. Phys.: Condens. Matter, 2017, 29, 273002.
45. T. R. Jow, J. Allen, M. Marx, K. Nechev, B. Deveney and S. Rickman, ECS Trans., 2010, 25, 3–12.
46. K. Tasaki, A. Goldberg, J.-J. Lian, M. Walker, A. Timmons and S. J. Harris, J. Electrochem. Soc., 2009, 156, A1019.
47. M. T. Dunstan, J. M. Griffin, F. Blanc, M. Leskes and C. P. Grey, J. Phys. Chem. C, 2015, 119, 24255–24264.
48. J. Christensen and J. Newman, J. Electrochem. Soc., 2004, 151, A1977.
49. Y. C. Chen, C. Y. Ouyang, L. J. Song and Z. L. Sun, J. Phys. Chem. C, 2011, 115, 7044–7049.
50. K. Xu, A. Von Cresce and U. Lee, Langmuir, 2010, 26, 11538–11543.
51. J. Mizusaki, H. Tagawa, K. Saito, K. Uchida and M. Tezuka, Solid State Ionics, 1992, 53–56, 791–797.
52. M. Dissanayake, Solid State Ionics, 1986, 21, 279–285.
53. N. Hornsveld, B. Put, W. M. M. Kessels, P. M. Vereecken and M. Creatore, RSC Adv., 2017, 7, 41359–41368.
54. J. Pan, Q. Zhang, X. Xiao, Y.-T. Cheng and Y. Qi, ACS Appl. Mater. Interfaces, 2016, 8, 5687–5693.
55. S. Shi, Y. Qi, H. Li and L. G. Hector, J. Phys. Chem. C, 2013, 117, 8579–8593.
56. J. Christensen and J. Newman, Proc. - Electrochem. Soc., 2003, 20, 85–94.
57. A. K. Rappe and W. A. Goddard, J. Phys. Chem., 1991, 95, 3358–3363.
58. Z. Liu, Y. Qi, Y. X. Lin, L. Chen, P. Lu and L. Q. Chen, J. Electrochem. Soc., 2016, 163, A592–A598.
59. A. Hajibabaei, C. W. Myung and K. S. Kim, Phys. Rev. B: Condens. Matter Mater. Phys., 2021, 103, 1–7.
60. A. Hajibabaei and K. S. Kim, J. Phys. Chem. Lett., 2021, 12, 8115–8120.
61. D. K. J. Lee, T. L. Tan and M.-F. Ng, ACS Appl. Mater. Interfaces, 2024, 16, 64364–64376.

---