Kinetic Monte Carlo Simulations Research Articles

Understanding Substrate Effects on 2D MoS2 Growth: A Kinetic Monte Carlo Approach

AbstractControlling the morphology of 2D transition metal dichalcogenides (TMDs) plays a key role in their applications. Although chemical vapor deposition can achieve wafer‐scale growth of 2D TMDs, a comprehensive theoretical framework for effective growth optimization is lacking. Atomistic modeling methods offer a promising approach to delve into the intricate dynamics underlying the growth. In this study, kinetic Monte Carlo (kMC) simulations are employed to identify crucial parameters that govern the morphology of MoS2 flakes grown on diverse substrates. The simulations reveal that large adsorption rates significantly enhance growth speed, which however necessitates rapid edge migration to achieve compact triangles. Substrate etching can tune the adsorption–desorption process of adatoms and enable preferential growth within a specific substrate region, controlling the flake morphology. This study unravels the complex dynamics governing 2D TMD morphology, offering a theoretical framework for decision‐making in the design and optimization of TMD synthesis processes.

Modeling of tungsten impurity transport and distribution in EAST based on multi-fluid and kinetic Monte Carlo simulations

In recent years, the materials of plasma facing components, such as divertor target plates, domes, and outer walls of tokamaks, such as ASDEX Upgrade, WEST, JET, EAST, and ITER, have been changed from carbon to tungsten because of its lower erosion and tritium retention rates. Impurities are produced by interactions between the plasma and the first wall. This study provides an investigation to simulate the transport and distribution of tungsten impurities in the edge plasma on EAST. The 2D multi-fluid edge plasma transport code SOLPS-ITER and 2D kinetic Monte Carlo impurity transport code DIVIMP were used in the simulations. The multi-fluid model in SOLPS-ITER and the kinetic Monte Carlo model in DIVIMP were employed to treat tungsten impurity ions. The 2D density contour distributions in the computational region and the 1D density radial profiles at the inner and outer midplanes of tungsten impurity particles with ionization states (W0–W+74) and the total tungsten particles with all charge states were obtained. With the heating power 1.5 MW and the line-averaged plasma density 2 × 1019 m−3, the results from SOLPS-ITER and DIVIMP show that the maximum density of tungsten ion with single ionization state is about 1014 m−3 and the total density of tungsten impurities with all charge states is about 1015 m−3 at the core boundary. To the best of our knowledge, the simulation results from SOLPS-ITER and DIVIMP are compared for the first time to benchmark SOLPS-ITER with the multi-fluid mode and DIVIMP with the kinetic model for tungsten impurity transport. The density distributions of tungsten impurities with different ionization states from SOLPS-ITER and DIVIMP are highly similar, and good agreement can be found under similar conditions involved in the calculation. From the comparison benchmark between SOLPS-ITER and DIVIMP for tungsten impurity transport, it can be concluded that the impurity transport approximation used by DIVIMP is good.

Diffusion-Limited Kinetics in Reactive Systems.

A proper representation of chemical kinetics is vital to understanding, modeling, and optimizing many important chemical processes. In liquid and surface phases, where diffusion is slow, the rate at which the reactants diffuse together limits the overall rate of many elementary reactions. Commonly, the textbook Smoluchowski theory is utilized to estimate effective rate coefficients in the liquid phase. On surfaces, modelers commonly resort to much more complex and expensive Kinetic Monte Carlo (KMC) simulations. Here, we extend the Smoluchowski model to allow the diffusing species to undergo chemical reactions and derive analytical formulas for the diffusion-limited rate coefficients for 3D, 2D, and 2D/3D interface cases. With these equations, we are able to demonstrate that when species react faster than they diffuse they can react orders of magnitude faster than predicted by Smoluchowski theory, through what we term "the reactive transport effect". We validate the derived steady-state equations against particle Monte Carlo (PMC) simulations, KMC simulations, and non-steady-state solutions. Furthermore, using PMC and KMC simulations, we propose corrections that agree with all limits and the computed data for the 2D and 2D/3D interface steady-state equations, accounting for unique limitations in the associated derived equations. Additionally, we derive equations to handle couplings between diffusion-limited rate coefficients in reaction networks. We believe these equations should make it possible to run much more accurate mean-field simulations of liquids, surfaces, and liquid-surface interfaces accounting for diffusion limitations and the reactive transport effect.

Precipitation kinetics of Cu in FeCu and FeCuX (X=B, N, and BN) alloys using kinetic Monte Carlo simulations

We employed kinetic Monte Carlo simulations to investigate the precipitation of Cu in FeCu and FeCuX (X=B, N, and BN) alloys. The parameters of kinetic Monte Carlo simulations were obtained through first-principles calculations. In this study, we utilized a vacancy-mediated model that considers the migration barriers dependent on local impurity atomic environments. The dynamical stability of both FeCu and FeCuX (X=B, N, and BN) alloys was examined by the virtual crystal approximation method. Through the assessment of advancement factor, cluster number density, and mean cluster radius, the simulation results indicated that both B and N promote precipitation of Cu. While B doping has a minor effect, N doping significantly promotes Cu precipitation.

Unified design flow for facilitating fast Li kinetics in layered oxide cathodes

Fast charging is desirable for the global use of full electric vehicles (EVs); however, it faces challenges based on Li kinetics. A structural design based on the electrostatic repulsion (ER) between cations over the full charging process is proposed herein to facilitate the fast Li kinetics of LiMO2 (M = transition metal) layered oxides in lithium-ion batteries (LIBs). This process is divided into two steps under the Li kinetic mechanism, in which LiO4 is a transition state upon migration. First, changing the O3- to O2-type layered cathodes causes the expansion of the LiO2 slab and elimination of the ER between Li and M ions in the hidden space based on the edge-shared mode for the Li-abundant region. Second, expanding the LiO4 space to release the stressed environment provides easy accommodation of Li ions through the ionic radius of the dopants in the Li-deficient region. Using unified atomistic calculations of density functional theory and kinetic Monte Carlo simulations, the systematic design flow can be understood in terms of the structural parameters, site energetics for the kinetic barrier, and statistical transitions based on various Co-based LiMO2 models. Our concrete understanding of local-to-macroscopic structural scales suggests a generalized direction for facilitating fast Li kinetics for realizing fast-charging LIBs.

Thermodynamics and kinetics of interconversion between platinum nanoparticles and cations in zeolites

In metal-containing zeolites, sintering and redispersion processes exercise control over the identity of metal site structures, but the thermodynamic and kinetic factors that influence these molecular-level processes are not completely understood. Here, we assess the ability of first principles informed free energy models (for supported and unsupported nanoparticles) and kinetic models integrating Ostwald ripening with atom trapping to describe the interconversion between Pt cations and nanoparticles encapsulated in chabazite (CHA) zeolites. Density functional theory-derived thermodynamic phase diagrams show that the interconversion between cations, favored in oxidizing environments, and particles, favored in reducing environments, is fully reversible within a wide range of their respective conditions (temperatures and pressures) for CHA and several other zeolite topologies. Kinetic Monte Carlo simulations of Pt redispersion are consistent with experimentally observed redispersion kinetics of encapsulated Pt nanoparticles in CHA zeolites, and model results suggest the zeolite host imparts additional stability for Pt nanoparticles. We envision our thermodynamic and kinetic models for Pt-CHA are also capable of describing nanoparticle and cation interconversions for other zeolite frameworks under similar conditions.

Kinetic coevolutionary models predict the temporal emergence of HIV-1 resistance mutations under drug selection pressure

Drug resistance in HIV type 1 (HIV-1) is a pervasive problem that affects the lives of millions of people worldwide. Although records of drug-resistant mutations (DRMs) have been extensively tabulated within public repositories, our understanding of the evolutionary kinetics of DRMs and how they evolve together remains limited. Epistasis, the interaction between a DRM and other residues in HIV-1 protein sequences, is key to the temporal evolution of drug resistance. We use a Potts sequence-covariation statistical-energy model of HIV-1 protein fitness under drug selection pressure, which captures epistatic interactions between all positions, combined with kinetic Monte-Carlo simulations of sequence evolutionary trajectories, to explore the acquisition of DRMs as they arise in an ensemble of drug-naive patient protein sequences. We follow the time course of 52 DRMs in the enzymes protease, RT, and integrase, the primary targets of antiretroviral therapy. The rates at which DRMs emerge are highly correlated with their observed acquisition rates reported in the literature when drug pressure is applied. This result highlights the central role of epistasis in determining the kinetics governing DRM emergence. Whereas rapidly acquired DRMs begin to accumulate as soon as drug pressure is applied, slowly acquired DRMs are contingent on accessory mutations that appear only after prolonged drug pressure. We provide a foundation for using computational methods to determine the temporal evolution of drug resistance using Potts statistical potentials, which can be used to gain mechanistic insights into drug resistance pathways in HIV-1 and other infectious agents.

Reverse Engineering of Radical Polymerizations by Multi-Objective Optimization.

Reverse engineering is applied to identify optimum polymerization conditions for the synthesis of polymers with pre-defined properties. The proposed approach uses multi-objective optimization (MOO) and provides multiple candidate polymerization procedures to achieve the targeted polymer property. The objectives for optimization include the maximal similarity of molar mass distributions (MMDs) compared to the target MMDs, a minimal reaction time, and maximal monomer conversion. The method is tested for vinyl acetate radical polymerizations and can be adopted to other monomers. The data for the optimization procedure are generated by an in-house-developed kinetic Monte-Carlo (kMC) simulator for a selected recipe search space. The proposed reverse engineering algorithm comprises several steps: kMC simulations for the selected recipe search space to derive initial data, performing MOO for a targeted MMD, and the identification of the Pareto optimal space. The last step uses a weighted sum optimization function to calculate the weighted score of each candidate polymerization condition. To decrease the execution time, clustering of the search space based on MMDs is applied. The performance of the proposed approach is tested for various target MMDs. The suggested MOO-based reverse engineering provides multiple recipe candidates depending on competing objectives.

Quantifying the energy of vacancy mediated diffusion of solute Cu, Mg, Ag and Zn atoms in an aluminum Matrix: An insight from ab initio calculations

In this work, a consistent and accurate set of energy barrier values for the migration process, assisted by a vacancy, of atoms of different chemical elements belonging to the two most relevant families of age-hardenable Al-Cu and Al-Zn-based alloys is given. These values were calculated using the Nudged Elastic Band method within the Density Functional Theory framework for a vacancy diffusing in simple configurations of solute atoms, specifically the main alloying elements Cu and Zn involving a single solute atom or more complex configurations, such as pairs of elements of the same or different chemical species. In both cases, in our calculations, the presence of Mg and/or Ag solute atoms is also considered. The results obtained are analyzed in terms of the energetic interaction among vacancies and the different solute atoms. The advantage of having a set of vacancy diffusion energy barrier values calculated under the same criteria is that it can be used as reliable input parameters for future kinetic Monte Carlo simulations to study the precipitation kinetics in a wide type of age-hardenable alloys.

Multiscale Investigation of the Mechanism and Selectivity of CO2 Hydrogenation over Rh(111).

CO2 hydrogenation over Rh catalysts comprises multiple reaction pathways, presenting a wide range of possible intermediates and end products, with selectivity toward either CO or methane being of particular interest. We investigate in detail the reaction mechanism of CO2 hydrogenation to the single-carbon (C1) products on the Rh(111) facet by performing periodic density functional theory (DFT) calculations and kinetic Monte Carlo (kMC) simulations, which account for the adsorbate interactions through a cluster expansion approach. We observe that Rh readily facilitates the dissociation of hydrogen, thus contributing to the subsequent hydrogenation processes. The reverse water-gas shift (RWGS) reaction occurs via three different reaction pathways, with CO hydrogenation to the COH intermediate being a key step for CO2 methanation. The effects of temperature, pressure, and the composition ratio of the gas reactant feed are considered. Temperature plays a pivotal role in determining the surface coverage and adsorbate composition, with competitive adsorption between CO and H species influencing the product distribution. The observed adlayer configurations indicate that the adsorbed CO species are separated by adsorbed H atoms, with a high ratio of H to CO coverage on the Rh(111) surface being essential to promote CO2 methanation.

Lignin structure dynamics: Advanced real-time molecular sensing strategies

Lignin, with its abundant and high energy density, offers potential for sustainable biofuels and chemicals. However, current research faces two primary challenges: limited experimental sampling frequency, impeding the real-time analysis and a lack of theoretical activation energies to elucidate thermodynamic and kinetic properties associated with structural changes under various processing conditions including reaction temperature (343 to 363 K), time (15 to 30 min), and chip size (1.0 to 5.0 mm). To address these challenges, we introduce the “real-time molecular-level sensing approach,” which considers delignification, redeposition, and de/repolymerization over time and various scales, leveraging sparse experimental data and theoretical activation energies. Our approach combines density functional theory, ab-initio molecular dynamics, and kinetic Monte Carlo simulations to deeply explore dynamic nature of lignin. Specifically, we gain insights into lignin content within biomass, molecular weight distribution, and monolignol ratio, improving lignin’s utilization during processing. This work establishes a new standard for comprehensive lignin analysis, fostering efficient lignin utilization by tracking their properties at the molecular level in the real-time manner.

Improved Approach for ab Initio Calculations of Rate Coefficients for Secondary Reactions in Acrylate Free-Radical Polymerization.

Secondary reactions in radical polymerization pose a challenge when creating kinetic models for predicting polymer structures. Despite the high impact of these reactions in the polymer structure, their effects are difficult to isolate and measure to produce kinetic data. To this end, we used solvation-corrected M06-2X/6-311+G(d,p) ab initio calculations to predict a complete and consistent data set of intrinsic rate coefficients of the secondary reactions in acrylate radical polymerization, including backbiting, β-scission, radical migration, macromonomer propagation, mid-chain radical propagation, chain transfer to monomer and chain transfer to polymer. Two new approaches towards computationally predicting rate coefficients for secondary reactions are proposed: (i) explicit accounting for all possible enantiomers for reactions involving optically active centers; (ii) imposing reduced flexibility if the reaction center is in the middle of the polymer chain. The accuracy and reliability of the ab initio predictions were benchmarked against experimental data via kinetic Monte Carlo simulations under three sufficiently different experimental conditions: a high-frequency modulated polymerization process in the transient regime, a low-frequency modulated process in the sliding regime at both low and high temperatures and a degradation process in the absence of free monomers. The complete and consistent ab initio data set compiled in this work predicts a good agreement when benchmarked via kMC simulations against experimental data, which is a technique never used before for computational chemistry. The simulation results show that these two newly proposed approaches are promising for bridging the gap between experimental and computational chemistry methods in polymer reaction engineering.

Island density scaling in reversible submonolayer growth with attachment barriers.

Island nucleation and growth play an important role in thin-film growth. One quantity of particular interest is the exponent χ, which describes the dependence of the saturation island density N_{sat}∼(D_{h}/F)^{-χ} on the ratio D_{h}/F of the monomer hopping rate D_{h} to the deposition rate F. While standard rate equation(RE) theory predicts that χ=i/(i+2) (where i is the critical island size), more recently it has been predicted that in the presence of a strong barrier to the attachment of monomers to islands, a significantly larger value χ=2i/(i+3) may be observed. While this prediction has recently been tested using kinetic Monte Carlo simulations for the case of irreversible growth corresponding to i=1, it has not been tested for the case of reversible island growth corresponding to i>1. Here we present a mean-field self-consistent RE method which we have used to study the dependence of the effective value of χ on D_{h}/F and barrier-strength for i=1,2,3, and 6. Both the no-nucleation-barrier case in which there exists a barrier for monomers to attach to islands larger than the critical island size (but not to smaller islands) and the nucleation-barrier case in which there is a barrier for monomers to attach to islands of all sizes are studied. In all cases, we find that the existence of attachment barriers significantly increases the effective value of χ for a given barrier strength. In addition, for i=1 we find good agreement between our extrapolated asymptotic value of χ and the theoretical strong-barrier prediction both with and without a nucleation barrier. In contrast, for i>1 the value of χ is significantly larger in the presence of a nucleation barrier than in its absence. In particular, while an asymptotic analysis of our results for i>1 also leads to excellent agreement with the strong barrier prediction in the presence of a nucleation barrier, in the absence of a nucleation barrier the asymptotic values are significantly lower.

Predicting pressure-dependent rate constants for the furan + OH reactions and their impact under tropospheric conditions.

In this work, we have evaluated the influence of temperature and pressure on the mechanism of furan oxidation by the OH radical. The stationary points on the potential energy surface were described at the M06-2X/aug-cc-pVTZ level of theory. In the kinetic treatment at the high-pressure limit (HPL), we have combined the multistructural canonical variational theory with multidimensional small-curvature tunneling corrections and long-range transition state theory. The system-specific quantum Rice-Ramsperger-Kassel theory was employed to estimate the pressure-dependent rate. In the HPL, the OH addition on the α carbon is the dominant pathway in the mechanism, producing a product via the ring-opening process, also confirmed by the product branching ratio calculations. The overall rate constant, obtained by a kinetic Monte Carlo simulation, reads the form koverall=5.22×10-13T/3001.10⁡exp1247(K/T) and indicates that the furan oxidation by OH radicals is a pressure-independent reaction under tropospheric conditions.

Two- and Three-Dimensional Recrystallization of Discrete Amorphous in C3H5-Molecular-Ion-Implanted Silicon Surface Analyzed by TCAD Simulation

Technology computer-aided design (TCAD) kinetic Monte Carlo simulations revealed the unique recrystallization processes of discrete amorphous regions connected to a buried amorphous layer in a C3H5-molecular-ion-implanted silicon (Si) substrate. The faithful simulation models show that the discrete amorphous regions are first recrystallized two-dimensionally in the lateral direction from both sides and separated from the buried amorphous layer. Then, the separated discrete amorphous regions are recrystallized three-dimensionally in the lateral and vertical directions from both sides and the bottom. We found that the first two-dimensional recrystallization of discrete amorphous regions is caused by the retardation of solid-phase epitaxial growth at the Si substrate surface and near the buried amorphous layer. We also found that the large (small) discrete amorphous regions require a long (short) two-dimensional recrystallization before separating from the buried amorphous layer. The transition point in the recrystallization dimension can be determined from the lateral recrystallization length and the equivalent radius of discrete amorphous regions.

Tailoring Hierarchical Structure and Rare Earth Affinity of Compositionally Identical Polymers via Sequence Control.

Macromolecule sequence, structure, and function are inherently intertwined. While well-established relationships exist in proteins, they are more challenging to define for synthetic polymer nanoparticles due to their molecular weight, sequence, and conformational dispersities. To explore the impact of sequence on nanoparticle structure, we synthesized a set of 16 compositionally identical, sequence-controlled polymers with distinct monomer patterning of dimethyl acrylamide and a bioinspired, structure-driving di(phenylalanine) acrylamide (FF). Sequence control was achieved through multiblock polymerizations, yielding unique ensembles of polymer sequences which were simulated by kinetic Monte Carlo simulations. Systematic analysis of the global (tertiary- and quaternary-like) structure in this amphiphilic copolymer series revealed the effect of multiple sequence descriptors: the number of domains, the hydropathy of terminal domains, and the patchiness (density) of FF within a domain, each of which impacted both chain collapse and the distribution of single- and multichain assemblies. Furthermore, both the conformational freedom of chain segments and local-scale, β-sheet-like interactions were sensitive to the patchiness of FF. To connect sequence, structure, and target function, we evaluated an additional series of nine sequence-controlled copolymers as sequestrants for rare earth elements (REEs) by incorporating a functional acrylic acid monomer into select polymer scaffolds. We identified key sequence variables that influence the binding affinity, capacity, and selectivity of the polymers for REEs. Collectively, these results highlight the potential of and boundaries of sequence control via multiblock polymerizations to drive primary sequence ensembles hierarchical structures, and ultimately the functionality of compositionally identical polymeric materials.

Unraveling the role of single-atom metal dopant over CeO2 catalyst in tuning catalytic performance of ethane oxidative dehydrogenation

Single-atom metal doped CeO2 (M1@CeO2) catalysts have exhibited better catalytic performance in methane and propane dehydrogenation, which may be applied to ethane dehydrogenation, revealing structure-activity relationship of M1@CeO2 catalysts is of great importance. In this work, DFT + U calculations together with kinetic Monte Carlo simulations were implemented to investigate the oxidative dehydrogenation stage (ODH) in ethane chemical looping-oxidative dehydrogenation (CL-ODH) over thirteen types of M1@CeO2 catalysts (M = Cr, Mn, Fe, Co, Ni, Ga, Mo, Ru, Rh, Pd, In, Ir and Pt). The results show that the doped single-atom metal and its adjacent O atoms act as active sites over four catalysts Rh, Pd, Ir and Pt doped CeO2, while only O atoms act as active sites over other nine catalysts M1@CeO2. Rh1@CeO2 was screened out to perform higher C2H4(g) selectivity of 98.8 % and C2H4(g) formation activity of 8.72×103 s−1 site−1 at 823.15 K with C2H6(g) partial pressure <0.2 bar, which is superior to the previously reported catalysts Pt, Pt3Sn, Pt3Ir, Ni3Mo, Co/g-C3N4 and Rh/ZrO2, and so on. A quantitative and quick descriptor, the binding strength of C2H5* intermediate is proposed to evaluate C2H4(g) formation activity. This study unravels the role of doped single-atom metal over CeO2 catalyst and help rationally design high-efficiency alkane dehydrogenation.

Adsorption and evolution of hydrogen molecules on hexagonal boron nitride monolayer: a combined DFT and kinetic monte-carlo simulations study

Density functional theory (DFT) and kinetic Monte-Carlo (kMC) simulation code has been combinedly used to study the adsorption and evolution dynamics of hydrogen molecules over a hexagonal boron nitride (h-BN) monolayer. Maximum adsorption energy from van der Waals curve is predicted to be around 60 to 70 meV using two different DFT functionals. Repulsive lateral interaction between two hydrogen molecules plays a key role in determining the maximum number of adsorptions inside one unit cell of h-BN. Bader charge analysis, electron localization function (ELF), total and partial density of states (DOS) plots have been included to understand the weak interaction going on between the adsorbent and substrate. The input energy parameters from the DFT calculation has been used to perform the kMC simulation for describing the adsorption, desorption and the diffusion pattern of hydrogen molecules with a given time of exposure to an empty h-BN substrate along with the overall surface coverage.

A DFT-based kinetic Monte Carlo simulation of multiphase oxide-metal thin film growth

Functional thin films of nanoscale metal pillars in oxide or nitride matrices known as vertically aligned nanocomposite (VAN) have gained much interest owing to their unique strain-coupled and highly anisotropic properties. So far, the deposition of these films has been explored mostly experimentally. In this work, a density functional theory (DFT)-based kinetic Monte Carlo simulation model using Bortz–Kalos–Lebowitz algorithm was developed to understand the growth of VAN films deposited by pulsed laser technique on mismatching substrates. The model has been parameterized and applied to understand the kinetics of growth thin films consisting of Au pillars in CeO2 matrix deposited on SrTiO3 substrates. The effects of pulsed laser deposition (PLD) conditions including the pulse frequency, deposition flux, and substrate temperature were explored. The simulations indicate that the Au pillar size and shape exhibit significant dependence on the PLD conditions. Namely, increasing the temperature increases the average pillar size and lowers the pillar density, and vice versa. In addition, the simulations revealed that increasing the deposition rate results in lowering the average pillar size and increasing the density. Particularly, the DFT results suggest that Au pillar size can be tuned during the initial growth of the first monolayer due to the significantly low activation barrier. Our analysis showed that the relationship between the average pillar size and pillar density is influenced by the kinetics. Furthermore, autocorrelation analysis showed that pillars self-organize in quasi-ordered patterns at certain windows of the deposition conditions, which is attributed to the complex nature of the chemical interactions in the system, the kinetics, and the deposition parameters.

Study of lithium transport in Li2O component of the solid electrolyte interphase in lithium-ion batteries

In lithium-ion batteries (LIBs), as a promising energy storage device, materials with fast lithium (Li) transport are required for high-power applications such as electric vehicles. Thus, a deeper understanding of Li transport in components of LIBs is crucial for improving their rate capability. In this study, the Li transport in lithium oxide (Li2O), as one of the key components of the solid electrolyte interphase (SEI) layer, was examined by a multiscale computational approach ranging from density functional theory (DFT) to Monte Carlo simulations. The DFT calculations were used to investigate the recombination of Frenkel pairs, their first-principles total energies, and the Li diffusion mechanisms. The effect of atomic configurations on both first-principles total energies and diffusion barrier energies was formulated by periodic and local cluster expansions. These formalisms were then incorporated into the Monte Carlo and Kinetic Monte Carlo (KMC) simulations to calculate the diffusion coefficient of Li. Our calculations revealed that the vacancy-mediated jump along the 〈100〉 direction in the antifluorite structure of Li2O possesses the lowest barrier energy compared to other diffusion mechanisms. The KMC simulations indicated that the diffusion coefficient of Li better converged with the direct experimental measurement when the recombination of Frenkel pairs was integrated into the simulations. At a temperature of 300 K, the KMC simulation yielded a Li diffusion coefficient of 3.8×10-12cm2/s in Li2O. This is only one order of magnitude larger than indirect experimental measurement, suggesting the accuracy of our formalism. Thus, our formalism for studying Li transport in Li2O will pave the path to a rational design of inorganic SEI in the future development of LIBs.