Kinetically corrected Monte Carlo-molecular dynamics simulations of solid electrolyte interphase growth

We present a kinetic approach to the Monte Carlo-molecular dynamics (MC-MD) method for simulating reactive liquids using nonreactive forcefields. A graphical reaction representation allows definition of reactions of arbitrary complexity, including their local solvation environment. Reaction probabilities and molecular dynamics (MD) simulation times are derived from ab initio calculations. Detailed validation is followed by studying the development of the solid electrolyte interphase (SEI) in lithium-ion batteries. We reproduce the experimentally observed two-layered structure on graphite, with an inorganic layer close to the anode and an outer organic layer. This structure develops via a near-shore aggregation mechanism. Introduction Understanding and controlling the chemical degradation of liquid electrolytes is key to enhancing the durability of lithium-ion batteries (LIBs). Even inside its electrochemical stability window, the electrolyte can reduce near the anode during charging. This initiates a cascade of reactions and culminates in the irreversible growth of a passivation layer, the solid electrolyte interphase (SEI). This layer is composed of various inorganic and organic species that contain lithium. The resulting capacity loss can be 10% in the first charge-discharge cycle and increases throughout the battery’s lifetime. The SEI is electronically insulating but ionically conducting and protects the electrolyte from further degradation while allowing lithium transport. The properties of the SEI and the overall battery performance depend on the initial electrolyte composition and its reactions near the anode. Successful simulations of SEI growth must provide a link between the molecular processes and the emergence of a microscopic structure. In this context, they can play a vital role in the in silico design of electrolytes for use in next-generation high performance LIBs. The general MC-MD hybrid method initially developed by Takenaka et al., later extended as the REACTER algorithm and independently by Biedermann et al., aims to build up microscopic structure by performing a suc-


Introduction
Understanding and controlling the chemical degradation of liquid electrolytes is key to enhancing the durability of lithium-ion batteries (LIBs). Even inside its electrochemical stability window, the electrolyte can reduce near the anode during charging. This initiates a cascade of reactions and culminates in the irreversible growth of a passivation layer, the solid electrolyte interphase (SEI). 1 This layer is composed of various inorganic and organic species 2 that contain lithium. The resulting capacity loss can be 10% in the first charge-discharge cycle and increases throughout the battery's lifetime. 3 The SEI is electronically insulating but ionically conducting and protects the electrolyte from further degradation while allowing lithium transport. The properties of the SEI and the overall battery performance depend on the initial electrolyte composition and its reactions near the anode. Successful simulations of SEI growth must provide a link between the molecular processes and the emergence of a microscopic structure. In this context, they can play a vital role in the in silico design of electrolytes for use in next-generation high performance LIBs.
The general MC-MD hybrid method initially developed by Takenaka et al.,[4][5][6][7] later extended as the REACTER algorithm 8-10 and independently by Biedermann et al.,11,12 aims to build up microscopic structure by performing a suc-cession of reactive cycles. These consist of reaction candidate searches and sampling, reaction execution and post-reaction MD components.
Candidates are identified and reactions executed using a pre-defined list of reactions. Separation of the reactive and dynamical steps in each cycle means that classical non-reactive MD can be used. Control over the sampling of candidates and length of MD being run is maintained, and reactions are guaranteed to be chemically viable. Similar methods have also been used to model the curing of composites 13 and could, in principle, be applicable to any reactive set of molecules.
Previous studies simulating SEI growth use a thermodynamic approach to MC-MD, accepting or rejecting reactive steps according to the Metropolis Monte Carlo scheme. Based on reaction energies, this ultimately leads to a Boltzmann distribution of reaction products. 14 Candidate sampling is often performed purely at random, and constant MD lengths are used in each reactive cycle. [5][6][7] This makes it difficult to simulate the competition between individual reactions occurring on different timescales.
Here, we present a kinetic approach to the MC-MD method in the form of kinetically corrected reaction candidate sampling and variable length post-reaction MD. Template reaction paths are pre-characterized using automated transition state searches based on density functional theory (DFT), 15 using the FlexTS module in BIOVIA Materials Studio. [16][17][18][19] We include spectator Li + ions and use an implicit solvent environment 20 to simulate the overall electrostatics. Corrections are based on the resulting activation barriers, and are applied to the candidate sampling scheme and postreaction MD lengths using equations inspired by the kinetic Monte Carlo method used in many branches of computational chemistry. 21 Our implementation uses the graphical and intuitive representation of template reaction paths in BIOVIA Materials Studio. 22 Coordination criteria and R-groups can be built into a reactant system of arbitrary complexity. These template reaction paths are then used to efficiently identify reaction candidates in simulation cells containing tens of thousands of atoms.
Combined with kinetically corrected MC-MD, this provides the link between individual molecular processes characterized with DFT and the emergence of a microscopic SEI structure produced with classical MD.
Our approach provides mesoscale insights into the structural dynamics of SEI growth. We provide a generalized efficient method to simulate the initial stages of SEI growth, capture condensation of an outer organic layer and the beginnings of crystallization of an inner inorganic layer. Furthermore, our method is able to reduce the computational cost of simulations relative to previous MC-MD implementations, thus providing the basis for a generalized and scalable electrolyte design workflow.

Introducing reactions into classical molecular dynamics
In our model anode/electrolyte system shown in Fig. 1a), identification of reaction candidates requires the matching of atom, bond and critical distance criteria (the ABC s) to the reactant of the corresponding template reaction path. These are constructed graphically in BIOVIA Materials Studio.
A template reactant can contain any number of fragments and is represented by a graph network where every atom forms a node, connected by edges corresponding to either bonds or specified critical distances. The ability to include general spectator species in template reactions allows coordination and micro-solvation criteria to be defined. This is particularly important when simulating SEI growth in LIBs, as the local coordination of Li + has a significant effect on the reaction energies and activation barriers of electrolyte degradation processes. 23 Our reaction definition also extends to redox reactions and those involving R-groups. Fig. 1b) shows an example template reaction path. It involves the dimerization of openethylene carbonate (oEC •− ) into ethylene dicarbonate (EDC 2− ) and ethene gas (C 2 H 4 ), both of which are common SEI components. Re- Figure 1: (a) model anode-electrolyte system, with the 1 nm thick reduction zone indicated by the green band adjacent to the anode. (b) example template reaction for the dimerization of the radical anion oEC •− to form EDC 2− and ethene gas. Critical distance criteria between reacting and coordinating atoms are shown in dashed lines. Atoms are colored by element (carbon -gray, oxygen -red, hydrogen -white, lithium -purple).
action occurs between the carbon radical of one oEC •− and the ethereal oxygen (O e ) of the other, and each oEC •− coordinates to Li + via their respective carbonate carbons (C c ). These reacting and coordinating criteria are defined through explicit specification of critical distances (dashed lines).
When searching for reaction candidates in situ, a pair of reacting or coordinating atoms, A and B, are within the critical distance, r crit , if their separation is less than the sum of their van der Waals radii. This can be scaled by a factor f VDW to control number of possible neighboring atoms to consider: Once identified and selected, a reaction candidate can be virtually reacted in situ by modifying bond orders, atomic charges and forcefield types to match those of the template product. Redox reaction candidates are only considered if they satisfy additional criteria. Here the center-of-mass of each reactant molecule of a reduction reaction must be within our defined reduction zone, as indicated by the green band in Fig. 1a). In addition to the appropriate bond order changes, the reducing electron is treated implicitly as an additional negative charge. This is introduced upon updating atomic parameters to those of the template product. This process can also be applied to the removal of electrons in oxidation reactions. Using bond orders and the overall charges of the reactants and products defined in our templates, we automatically detect the charge transfer character of each reaction without requiring explicit user input.
In the present study, the thickness of the reduction zone has been set to r e − = 1 nm from the anode surface, which means that reduction reactions are restricted to that area of the simulation cell. Based on work by Lin et al. 24 that estimated the insulating thickness of a layer of Li 2 CO 3 as 3 nm (at which electron tunnelling probability is essentially zero 25 ), our value averages the distance at which appreciable electron tunnelling can occur.

A kinetic approach to Monte Carlo-molecular dynamics
The MC-MD method requires the execution of many hundreds of reactive cycles as outlined in Fig. 2. Our implementation allows both thermodynamic and kinetic sampling procedures using Metropolis Monte Carlo and kinetic Monte Carlo inspired approaches respectively. Prior to performing MC-MD, the model system is constructed using the Amorphous Cell module in Materials Studio and template reactions are defined (see Supporting Information).
A reactive cycle begins with a search of the model system for all reaction candidates corresponding to template reaction paths. If none are found, a short MD simulation of length τ update is run to update the total configuration and another search is performed (update loop in Fig. 2). Execution of 200 consecutive update loops results in a simulation being terminated. A new reactive cycle is initiated once at least one candidate is found, and a simulation is considered complete if sufficient reactive cycles have been performed.
Next we determine the selection probability for each candidate. Several studies using a thermodynamic approach have selected candidates uniformly at random. 5-7 The kinetic correction to the selection probability for each template path p is based on the reaction rate: Here k B is Boltzmann's constant, T = 298 K is the temperature of the system, A p = k B T /h is the pre-exponential factor, and E p a are the activation barriers obtained from DFT. A candidate is selected at random using weights based on Eqn. (2).
The total reaction rate is: and is based only on the current state of the system, e.g. the number of reaction candidates, N p cand , present for each template path and their associated rates. To avoid numerical instability in the post-reaction MD, local optimization of all atoms involved in the reaction is required.
Following the reaction step, an MD simulation of length τ evolves the system. The evolution time may be fixed throughout the simulation, τ const , as used in previous thermodynamic and kinetically corrected studies. 4-7,12 Here we vary this timescale according to the total reaction rate in Eq. (3) as: where R ∈ (0, 1] is another random number. Together with the selection process, our procedure corresponds to the kinetic Monte-Carlo approach. 21 As τ var can vary significantly according to the current state of the system, and many hundred reactive cycles are required, we limit τ var according to: with τ min = 100 fs and τ max = 20 ps. This effectively yields instantaneous reactions when barrierless reactions dominate the total rate, and significant diffusion when only highly activated processes are possible.

Electrolyte degradation reaction network
For reactive systems such as the anodeelectrolyte interface in LIBs, the complexity of the potential energy landscape allows many possible reactions at any given time. This requires simplification by creating a set of template reactions containing the most important processes. Construction of the network used here was context-driven based on previous studies. 23,26 The reaction paths for the degradation of EC (I -VIII) are shown in Fig. 3, along with the activation barriers calculated using DFT. All reactions are strongly exothermic (see Supporting Information). We have focused on EC due to its ubiquitous use in liquid electrolytes and active anode chemistry at typical voltages in battery cells. It is preferentially reduced over other carbonate solvents in widely used electrolyte due to its high dielectric constant and polarity. 1 Initially only the 1e reduction of EC (path I) can occur, which results in a change of hybridization of the carbonyl carbon (C c ) from sp 2 to a partial sp 3 character that stabilizes the newly formed radical. Predicted by Wang et al., and confirmed in our calculations, coordination of the carbonyl oxygen O c of EC to Li + aids reductive decomposition. 23 This reduction produces the closed radical anion cEC •− which is stabilized in a dielectric medium compared to the gas phase. It can undergo a ring opening (path II, E a = 9.85 kcal mol −1 ) via homolysis of the HC-O e (ethereal oxygen) bond to form the open radical anion oEC •− , where the unpaired electron is localized on the carbon of the terminal CH 2 group, and the negative charge delocalizes over the carbonate group at the other terminus. oEC •− is an intermediate for many termination reactions forming common SEI components. It can undergo a second reduction at the anode surface, fragmenting into ethene gas (C 2 H 4 ) and the inorganic dianion CO 2− 3 . It can also react with cEC •− via radical pairing and fission of the C c -O e bond in cEC •− to form the alkoxyester 23 (O(CH 2 ) 2 CO 2 (CH 2 ) 2 OCO 2 ) 2− with a relatively high barrier of 22.18 kcal mol −1 . oEC •− can also undergo two self-dimerization reactions. Barrierless pairing of radical carbon atoms (path V) forms butylene dicarbonate (BDC 2− ). Ethylene dicarbonate (EDC 2− ) can be formed formed either by the barrierless pairing of the radical carbon of one oEC •− and O e of another (path VII -shown as the example in Fig. 1), or the pairing of radical carbon and O c (path VI), which due to the significant electronic population on O c proceeds with a very high activation barrier of 40.5 kcal mol −1 .
Nucleophilic attack of CO 2− 3 on the electrophilic sp 3 carbon of EC (path VIII) forms EDC 2− . Further investigation shows that the Li + coordination state of reactant species of path VIII strongly influences the activation barrier. Coordination of EC to Li + reduces it (relative to uncoordinated EC) by providing a more effective electron sink that aids the concerted S N 2-like C-C bond formation and HC-O e fission, confirmed by population analysis.
Charge can delocalize onto Li + when it is coordinated with CO 2− 3 . This weakens the nucleophilicity of oxygen and thus increases the activation barrier. Due to the strong likelihood of CO 2− 3 being coordinated to at least one Li + cation in solution, we have chosen to require coordination of both EC and CO 2− 3 to Li + in path VIII, which proceeds with an activation barrier of 8.02 kcal mol −1 .

Larger r crit for increased competition and simulation lifespan
The computational cost of a simulation is driven by two factors: i) the number of searches required to find at least one reaction candidate, ii) the number of MD steps during the postreaction MD and/or the update loop, with the candidate sampling and selection, virtual reaction execution and local geometry optimization having a relatively minor cost-percycle. Both factors are linked. Any failure to find candidates results in execution of the update loop, adding more MD and another search. A larger distance r crit increases the likelihood of finding reaction candidates. Changing the van der Waals scale factor f VDW in Eqn. (1) therefore has a significant effect on simulation performance. Fig. 4 shows the total number of candidates found in each of the first 250 reactive cycles of a simulation, for 3 different starting configurations and different f VDW . Figure 4: The number of candidates found at each reactive cycle for 3 different van der Waals scaling factors; f VDW = 1.00 (red), 1.25 (blue), 1.50 (green). Simulations run for 3 different starting cells.
In the first reactive cycle, where each respective starting cell is identical in terms of atomic positions and bonding topology, a higher f VDW provides more candidates. This also influences sampling over time. Comparison of the curves shows that they do not simply scale with f VDW , but rather that they broaden as well. The number of candidates being found affects the competition in the sampling pool as well as the post-reaction MD lengths calculated with Eqn. (4), influencing the nanometre-scale structural dynamics of SEI growth over hundreds of reactive cycles (see Supporting Information). By increasing r crit , a more diverse candidate pool is maintained for longer, enhancing the ability of our stochastic Monte Carlo scheme to select for kinetically competitive reactions.
In addition, more reactive cycles can be executed whilst only requiring one candidate search per cycle, in turn reducing the reliance on the update loop. The overall computational cost of a simulation is therefore reduced by minimizing i) and ii).
A value of f VDW = 1.50 in Eqn.
(1) was found to be most appropriate for extending simulation lifespan, maintaining competition, and providing a more realistic picture of SEI growth. A value larger than this can result in unphysical bonds forming, difficulty in local structure relaxation and destruction of conformational information implicit in the specified critical distances. Fig. 5 shows a comparison between a typical simulation with τ const = 8 ps and one using variable MD time. The 8 ps simulations were chosen as they were most similar to the variable approach in terms of the number of candidates found, path selection probabilities and evolution of species counts (see Supporting Information for a full study of different τ const values). For the constant approach the maximum shared reactive cycle reached by the three best performing cells was ∼ 650. For the variable approach, this reactive cycle was ∼ 1200.

Variable length MD
In the first 650 reactive cycles, both approaches show similar trends in the number of candidates found (top row) and the selection probabilities (middle row), due to the use of kinetically corrected sampling. The resulting cascade of reactions and build-up of products Results are shown as averages over the three best performing (i.e. highest reactive cycle reached) of five different starting cells. Kinetically corrected selection was used with f VDW = 1.50. Top rowthe total number of candidates for all paths found at each reactive cycle (black line) and the rolling average of the number of candidate searches required per cycle (green line). Middle row -the path selection probabilities for the template reaction paths (I -VIII), block averaged every 10 reactive cycles and stacked. Final row -species counts as a function of reactive cycle. The primary axis measures the count of the initial reduction product, cEC •− , with the secondary axis for all other species. The standard deviation in species count across the three starting cells are indicated with error bars. Successive nanoseconds of total elapsed MD time are marked with vertical gray dashed lines.
(final row) is also similar. For the variable MD time, standard deviations in the species counts are smaller, particularly for cEC •− . The reaction cascade is less chaotic and more resilient to starting conditions due to the allocation of MD time. When the initial path selection probabilities are dominated by the barrierless path I, τ var calculated on-the-fly frequently falls below the τ min lower limit in Eqn. (5). These reductions are executed instantaneously in succession, with no MD being run between them.
The constant MD time approach cannot capture this effect, as it systematically overestimates the timescales of barrierless reactions. This unphysical allocation of MD time also increases the computational cost, as shown by inspecting the reactive cycle at which each successive nanosecond of total elapsed MD time occurs (gray dashed lines). At reactive cycle 315, where the variable approach reaches 1 ns of elapsed MD time, the constant approach has already elapsed ∼ 2.5 ns. By reactive cycle 650, the variable approach has elapsed 3.38 ns of MD time, 38% lower than the 5.47 ns elapsed by the constant approach.
The variable approach is better able to allocate MD time in the portion of the simulation where competition between candidates is high, better representing reactive timescales and doing so with reduced computational cost-percycle. It also accounts accurately for the diffusive dispersal of species between reactive steps. Naturally, as reactive events become rare, all simulations will cross over into a regime where competition is significantly reduced, candidates are harder to find and progress is driven with increasing reliance on the update loop. This can be seen by inspecting the rolling average of the number of candidate searches required per cycle (green line, top row in Fig. 5).
An increase of the average above 1 indicates a gradual progression into the rare-event regime, and is accompanied by a shortening of the intervals between each successive nanosecond of elapsed MD time. For the variable τ var approach, this transition is smoother and occurs at a higher reactive cycle compared to the constant approach. Factors i) and ii) introduced previously are minimized for longer, allowing for more extensive simulation of SEI growth within a given set of computational resources.
SEI growth via near-shore aggregation Fig. 6 shows SEI growth over the course of 1300 reactive cycles in our ∼ 24k atom graphite-EC/LiPF 6 model system. By reactive cycle 1300, the SEI consists of a dense ∼ 5 nm thick outer organic layer, as well as a thin inner inorganic layer deposited directly on the anode surface. This is in line with experimental and computational studies that reveal the SEI as a two-layered structure. 27,28 Beginning early in the simulation CO 2− 3 (red) forms a stable layer on the anode, bound by strong electrostatic forces to multi-coordination state Li + , which increases in mass density as the simulation progresses. This presents the start of crystallization of Li 2 CO 3 , identified by X-ray photoelectron spectroscopy (XPS) to be one of the main components of the inner inorganic layer. 29 Ethene gas (gray), one of most abundant gases trapped in the SEI regardless of the type(s) of carbonate the electrolyte is based on, 30 accumulates between the inorganic and organic layers.
Here we are able to track the structural dynamics of how the SEI forms. The first 200 reactive cycles are dominated by the initial reduction reaction (path I), which occurs mostly as instantaneous changes in bonding topology owing to the barrierless nature of the path. Except for a few cEC •− (yellow) that sparsely populate the anode surface, the production of a high density of negative charge leads to minor exfoliation of the bulk electrolyte from the surface.
Around cycle 300, aggregates of BDC 2− (blue) and EDC 2− (green) form, held together by networking Li + ions (purple) with a two-or three-fold coordination state. These diffuse relatively freely into the bulk electrolyte, away from the anode surface. As more dimerization reactions occur, larger aggregates form (see Supporting Information) and the coordination state of Li + increases to 4-or 5-fold. By reactive cycle 600 an organic layer starts to condense. Initially, the layer is sparsely populated and relatively close to the anode surface.
Neutral EC molecules are continuously able to approach the reduction zone throughout as they are less strongly bound in large organic aggregates linked by Li + and less repelled from the anode surface compared to larger and more polar species such as BDC 2− and EDC 2− . Once reduced, cEC •− diffuses away from the surface, is converted into oEC •− (path II) and dimerizes (paths V, VI, VII).
From reactive cycle 600 onwards, the density and thickness of the organic SEI layer increases. Simultaneously, we see an increase in the layer stability against dissolution or disintegration into the surrounding electrolyte, owing to the strongly bound matrix created by BDC 2− and EDC 2− coordinated to linking Li + (see Supporting Information). The layer also extends further into the bulk EC solution, with the greatest extent of the SEI from the anode surface found to be ∼ 13 nm. At the same time, the Li + mass density separates with little remaining strongly bound in the inorganic layer on the anode surface, but most having been consumed within the organic layer. PF − 6 is expelled into the bulk electrolyte.
The mechanism of growth by means of continuous EC reduction and organic aggregation in the region close to but away from the anode surface follows the 'near-shore aggregation mechanism' of SEI growth proposed by Ushirogata et al. 31 in their computational DFT-MD study. Contrary to the 'surface-growth' mechanism, this helps to explain SEI thicknesses spanning tens of nanometres, observed here and also experimentally. 32 Figure 6: SEI growth in our ∼ 24k atom graphite-EC/LiPF 6 model system. Unit cell structures and mass density profiles are shown every 100 reactive cycles. Reaction products are presented in CPK style (full vdW radius) and have been colored according to their identity, in accordance with the boxed colors in Fig. 3. In-tact EC and PF − 6 molecules are presented in line style for visual clarity. The mass density profiles of Li + have been scaled by a factor of 10. The distribution of the PF − 6 anions across the unit cell is visualized using the mass density profile of phosphorous.

Discussion
Our kinetic approach to MC-MD provides a generalized, physically correct, and computationally efficient method for simulating a twolayered SEI structure from a model system and list of template reactions. While not reproducing exact dynamics on the atomic scale, we capture the important structural dynamical features on a nanometre level. The outer organic layer of the SEI contains mostly BDC 2− and EDC 2− bound by Li + . It forms via the nearshore aggregation mechanism, where continuous reduction of intact EC molecules at the anode and dimerization reactions away from it produces a heterogeneous SEI structure spanning ∼ 13 nm. While our reaction network considers only EC degradation paths, we also observe the beginnings of Li 2 CO 3 crystallization on the anode surface as part of the inner inorganic layer. In-clusion of PF − 6 decomposition paths in the production of fluoride-containing compounds will provide a more complete picture of the growth of an insulating inorganic matrix, of which LiF is a major component. 29 Our reaction network is obtained by precharacterizing electrolyte degradation paths using automated transition-state searches based on DFT. The effects of explicit micro-solvation and implicit solvent can be captured and used to introduce kinetic corrections to our classical simulations through our generalized MC-MD method.
Our simulation scheme focuses on the physically correct timescales for both reactive and diffusive processes. In situations where multiple reactions compete, reaction paths with lower barriers are favoured over more activated paths. While species that can be produced are limited to those included in the reaction network, those that do form are ultimately dictated by activation barrier. When using uncorrected sampling, populations are more an artefact of the number of candidates found and appearances in the reaction network.
Variable length post-reaction MD involves onthe-fly calculation of reactive timescales, τ var , inspired by kinetic Monte Carlo. Depending on the number and activation barriers of the reaction candidates present in the system at a given instant, this approach can distinguish between instantaneous transitions and those activated over longer timescales. A more physical allocation of MD time is also computationally efficient, lowering the overall cost-per-cycle.
The critical distance, r crit , used to identify pairs of reacting species, can be controlled by scaling the van der Waals factor, f VDW . Increasing r crit extends the number of reactive cycles and lowers cost-per-cycle by reducing reliance on the update loop. It also enhances the sampling power of the kinetically corrected Monte Carlo selection scheme by maintaining competition for longer. While conformational information is implicitly encoded in the definition of critical distances, our tool can be extended to include bond angle and dihedral criteria.
Further optimization and parallelization of our candidate searching algorithm will allow us to simulate larger systems over longer timescales. With the main limitation of this method being the timescales that can be accessed by MD, more sophisticated biasing techniques could be used to enhance dynamical sampling of bonding topologies and treatment of simulations in the rare-event regime.

Conclusions
The atomistic reaction representation in BIOVIA Materials Studio allows graphical and intuitive construction of template reactions of arbitrary complexity. Micro-solvation criteria, redox character, and reactions involving Rgroups can also be defined. Using our generalized MC-MD approach with kinetic corrections, only the model system, list of template reaction paths and their activation barriers need to be defined to simulate SEI growth.
Going forward, the design and characterization of novel anode, electrolyte and additive systems can start with the exploration of reaction space using chemical reaction networks. 33,34 A list of the most important reactions likely to occur and contribute to SEI growth could be generated, with paths characterized using FlexTS and defined as template reaction paths using our reaction builder for input into the MC-MD method.
This provides the basis for an automated workflow, and could be applied to any reactive system. For battery performance, this could allow for screening of electrolyte compositions, with resulting SEI structures evaluated for their ability to prevent electrolyte degradation and avoid significant capacity loss. Combined with with workflows already implemented in Materials Studio for evaluating transport properties, 35 this enables the multi-objective optimization of electrolyte composition within the same computational framework. These optimizations, alongside consideration of other ageing processes such as lithium plating and electrode fracture, will play a key role in the development of next-generation high performance LIBs of greater sustainability and durability.

Computational Details
Model Systems were constructed as follows. The anode was modelled using a crystal of AB-stacked graphite, with layers lying in the xz plane and the surface in the xy plane at z = 0. A unit cell with (A, B, C) dimensions lying along the (x, y, z) axes respectively was created, containing the anode crystal and a vacuum region adjacent to it, extending along the C-dimension of the cell. The vacuum region was packed with a liquid electrolyte solution to a target density of 1.32 g cm −3 using Materials Studio Amorphous Cell. The electrolyte was comprised of ethylene carbonate (EC) solvent molecules and lithium hexafluorophosphate (Li + / PF − 6 ) ion pairs to a concentration of ∼ 1.0 mol dm −3 at which Li-ion conductivity in EC-based electrolyte is close to its maximum. 35 Geometry optimization and NVT-MD equilibration of the cell was performed using Materials Studio Forcite. Finally, the trivalent carbon atoms on the surface of the anode (at z = 0) were each assigned a charge so that the total surface charge per unit cell was -2 e, in order to create a potential difference of 2 V between anode and electrolyte. 5,6,36 The anode acts as a steric barrier between the electrolyte and its periodic images in the Cdirection, as well as a charged surface adjacent to which a reduction zone is defined. The anode was constrained at all points in construction and MC-MD simulations. Construction of model systems was repeated to generate configurationally distinct starting cells.
Results shown in Figs. 4 and 5 were generated for a unit cell of linear dimensions (34.08× 37.00× 124.60) A. The electrolyte region was 100 A in length (along the Cdimension) and contained 983 EC molecules and 74 Li + / PF − 6 ion pairs, corresponding to a 1 M salt concentration. The anode was comprised of 10 graphene layers of 3.70 A spacing (corresponding to lithiated graphite). The NVT-MD equilibration length was 30 ns.
Results in Fig. 6 were generated using the model system shown in Fig. 1a). The unit cell had linear dimensions (38.35 × 40.20 × 173.70) A. The electrolyte region was 150 A in length (along the C-dimension) and contained 1850 EC molecules and 140 Li + / PF − 6 ion pairs, again using a 1 M salt concentration. The anode was comprised of 12 graphene layers of 3.35 A spacing (corresponding to pure graphite). The reduced spacing was to prevent intercalation of species into the anode. NVT-MD equilibration length was 5 ns.
DFT Calculations of reaction energies and activation barriers were performed using DMol 3 15 in Materials Studio at the m-GGA/TPSS level of theory. Reaction barriers were calculated using the doublynudged elastic band (DNEB) method 17 as implemented in FlexTS/DMol 3 , with barrierless (E a = 0) reactions confirmed as such by observing a monotonically decreasing energy landscape along the reaction path. Calculations were run in a simulated solution-phase environment as implemented in COSMO, 20 using the extrapolated dielectric constant for pure ethylene carbonate at 298 K of = 95.3. 37,38 MD Simulations were run using Materials Studio Forcite 22 as a constant-number, -volume, and -temperature (NVT) ensemble, with an integration timestep of 1 fs and the temperature maintained at 298 K with a Berendsen thermostat. The COMPASSIII forcefield 39 was used, with PPPM electrostatic summation. Geometry optimizations were also performed with Materials Studio Forcite.
Acknowledgement The authors thank Andreas Heuer, Diddo Diddens and Victor Milman for helpful discussions. The van der Waals scaling factor, f VDW The van der Waals scaling factor, f VDW , controls the leniency of the critical distance criteria used when searching for reaction candidates. A higher f VDW value generally results in more reaction candidates being found (see Fig. 4 in the main text). Fig. S1 shows the species counts of simulations using different f VDW values (1.00, 1.25 and 1.50). In simulations using f VDW = 1.00 and 1.25, termination products of the reaction network form quickly as part of a surface growth mechanism, and removal of Li + from the reduction zone rapidly reduces the number of candidates being found. As a result, the simulations terminated early as 200 update loops (see Fig. 2 in the main text) are run consecutively without any reaction candidates being found. The use of a higher f VDW value affects the structural dynamics of SEI growth on the nanometre scale, influencing simulation performance and lifetime.

Supporting Information Available
Continuous reduction of electrolyte, and condensation of an organic layer away from the anode surface as part of a near-shore aggregation mechanism S1 produce longer lasting and S-1 better performing simulations. Growth of a heterogeneous SEI structure that spans tens of nanometres is observed, in line with experimental findings.

Kinetically corrected candidate sampling
Reaction candidates are sampled once identified using the stochastic Monte Carlo selection scheme described in the main text. The simplest way to do so is to select one at random, without applying any physical corrections (i.e. uniformly). Alternatively, each candidate's selection probability can be weighted using the molecular rate of the template reaction path is corresponds to prior to random selection. This weighting is done using pre-calculated activation barriers based on DFT. In this study, E a values shown in Fig. S3 were used.  produced early in the simulation. These organic components bind strongly to Li + and diffuse away from the reduction zone in aggregates. The reduction zone becomes stripped of Li + , making reduction of EC unlikely and candidates hard to find. As the termination products by definition do not react further, the uncorrected/constant simulations are terminated early.
Relative populations of termination products present at completion of uncorrected/constant simulations are an artefact of the number of template path the products appears in. EDC 2− , for example, produced by paths VI, VII and VIII, is dominant. BDC 2− on the other hand, produced only by the barrierless path V, is present in much smaller quantities. The alkoxyester is also present in amounts comparable to BDC 2− despite being produced only by the highly activated path III. were present in small but comparable quantities. The alkoxyester was not found in any quantity.
Uncorrected candidate searches on the other hand over-sample kinetically unfavourable reactions. This leads to unphysical accumulation of some species with uncontrolled effects on the overall SEI structure. Kinetically corrected sampling takes into account both the number of candidates found and their activation barriers. As a result, the initial stages of the simulation is dominated by reduction of EC. This provides a more physical treatment of the competition between paths in the anode-electrolyte system under kinetic control. In practice, we find that higher reactive cycles can be accessed while maintaining a high level of S-5 competition. The use of corrected selection reduces the overall computational cost-per-cycle by reducing reliance on the update loop and thereby minimizing factors i) and ii) introduced in the main text. It also contributes to a near-shore aggregation mechanism of SEI growth that produces a heterogeneous SEI structure spanning tens of nanometres.

Template reaction paths
Fig . S3 shows the graphical representation of the reactants and products of the eight template reaction paths considered in this study, as built using our reaction builder tool in BIOVIA Materials Studio. S12 Redox character, reaction energies and activation barriers are also indicated. All reactions where highly exothermic, hence reverse reactions were not considered, but the activation barriers varied significantly.

SEI Growth
Structures Unit cells for the structures shown in Fig. 6 in the main text are provided in Crystallographic Information File (CIF) format. Files are named according to the reactive cycle to which they correspond, i.e. 0.cif corresponds to the starting cell and 800.cif to the cell after reactive cycle 800 was executed.

Layer stability
In order to ascertain the relative stabilities of the various components of the SEI, cells at 100 reactive cycle intervals were extracted (i.e. the structures shown in Fig. 6 of the main text) and non-reactive molecular dynamics (of length 5 ns) were run on them. Fig. S4 shows the mass density profiles for the cells extracted after reactive cycles 100, 500, 900 and 1300, before (i.e. at 0 ns) and after (i.e. at 5 ns) the MD simulation, with the profile for S-6 Figure S3: The graphical representation of reactants and products for the eight template reaction paths considered in this study. Activation barriers, E a , and reaction energies, ∆E r , (both in kcal mol −1 ) for the paths are included, calculated using automated transition-state searches S7-S10 using DFT in BIOVIA Materials Studio DMol 3 . S11,S12 Atoms are colored by element (carbon -gray, oxygen -red, hydrogen -white, lithium -purple).
S-7 Figure S4: Mass density profile scatter matrix. Each column shows the mass density profiles of each species in the unit cell at the start (0 ns) and end (5 ns) of the non-reactive molecular dynamics simulation.
Only profiles for unit cells after reactive cycles 100, 500, 900 and 1300 are shown. Unit cell structures are extracted from simulations of the larger model system of length 15 nm, as shown in Fig. 6 in the main text.
The mass density profiles of Li + have been scaled by a factor of 10. The distribution of the PF − 6 anions across the unit cell is visualised using the mass density profile of phosphorous. each species shown in its own color. Each of these species profiles at 0 and 5 ns were then compared quantitatively by calculating their Pearson correlation coefficient. This provides a metric of similarity of the profiles along the Z-dimension of the unit cell (perpendicular to the anode surface), as shown in Fig. S5, indicating the stability of the various component layers to disintegration.

Consumption of Li +
In this study Li + acts only as a coordinating species; it is not chemically changed by any of the reactions. The proportion of Li + in the unit cell that has acted as a coordinating species in at least one reaction as a function of reactive cycle is shown in Figure S6b). As higher reactive cycles are reached and the condensation of separate inorganic and organic layers becomes more extensive, Li + ions become strongly bound within these layers. In the organic layer in particular, multi-coordination state lithium forms important bridges between organic components and thus is coordinated in place, available for participation in S-8 Figure S5: The Pearson correlation coefficient between the mass density profiles of each species at the start (i.e. 0 ns) and end (i.e. 5.0 ns) of a non-reactive dynamics simulation. These were run on the unit cell structures shown in Fig. 6 in the main text (i.e. every 100 reactive cycles). Figure S6: (a) A typical organic aggregate comprised of BDC 2− (blue), EDC 2− (green) and oEC •− (orange) fragments connected by multi-coordinate state Li + (purple). Extracted from the unit cell at reactive cycle 600, shown in Fig. 6 in the main text. (b) The proportion of Li + ions in the unit cell that have acted as a coordinating species at least once in the course of the production run simulation, as a function of reactive cycle. Data is extracted from simulations of the larger model system of length 15 nm. dimerization reactions. Large organic aggregates, such as the one shown in Fig. S6a), form as part of a near-shore aggregation mechanism of SEI growth. S1 S-9