---
author:
- Chiara Panosetti
- 'Simon B. Anniés'
- Cristina Grosu
- Stefan Seidlmayer
- Christoph Scheurer
bibliography:
- 'references.bib'
title: 'DFTB modelling of lithium intercalated graphite with machine-learned repulsive potential'
---

Abstract
========

Lithium ion batteries have been a central part of consumer electronics for decades. More recently, they have also become critical components in the quickly arising technological fields of electric mobility and intermittent renewable energy storage. However, many fundamental principles and mechanisms are not yet understood to a sufficient extent to fully realize the potential of the incorporated materials. The vast majority of concurrent lithium ion batteries make use of graphite anodes. Their working principle is based on intercalation—the embedding and ordering of (lithium-) ions in the two-dimensional spaces between the graphene sheets. This important process—it yields the upper bound to a battery’s charging speed and plays a decisive role for its longevity—is characterized by multiple phase transitions, ordered and disordered domains, as well as non-equilibrium phenomena, and therefore quite complex. In this work, we provide a simulation framework for the purpose of better understanding lithium intercalated graphite and its behaviour during use in a battery. In order to address the large systems sizes and long time scales required to investigate said effects, we identify the highly efficient, but semi-empirical Density Funtional Tight Binding (DFTB) as a suitable approach and combine particle swarm optimization (PSO) with the machine learning (ML) based Gaussian Process Regression (GPR) to obtain the necessary parameters. Using the resulting parametrization, we are able to reproduce experimental reference structures at a level of accuracy which is in no way inferior to much more costly *ab initio* methods. We finally present structural properties and diffusion barriers for some exemplary system states.

Introduction {#intro}
============

Within the past decade, studies investigating the consequences of man-made climate change [@Sharp2011; @Fisher2012; @Program2018] have become more specific, the predicted time frames shorter and the warnings more urgent. The immediate and radical reduction of carbon dioxide emissions by replacing fossil fuel based energy sources with renewable ones has been found to be the only reasonable approach to at least limit those consequences. [@Anderson2016] While the generation of electric energy from wind and sun is already quite advanced and efficient, its storage and transport are the main factors holding it back compared to coal and oil. Currently, two main approaches are being pursued in order to eliminate these drawbacks. One aims directly at the synthesis of alternative liquid or gas-phase fuels. The other intends to improve upon existing battery technology—especially lithium ion batteries—enough, to make it a serious contender in terms of energy sustenance. In this work, we intend to lay some groundwork for gaining deeper insight into some of the atomistic mechanisms limiting the (dis-)charging speed and lifetime of the most common types of lithium ion batteries, with graphite intercalation anodes. Ever since graphite was ascertained experimentally and theoretically to be an excellent candidate as an anode for Li-ion batteries, numerous attempts were made at fully describing the working system. [@Hennig1959; @Guerard1975; @Hawrylak1984; @Conard1994; @Nitta2015] Most of the electrochemical properties of the anode material itself are well-known. However, in particular transport processes during strongly driven operating conditions, like fast charging, are only poorly understood at a microscopic level. These technologically important macroscopic conditions are accompanied [*e.g.*]{} by temperature variations, leading to a capacity fade during ageing, as well as lithium plating. All of them limit the lifetime of the battery. [@Gallagher2016; @Wandt2018; @Yang2018] Against this background, experiments and theory are pushed quite far to gain insight into the real processes occuring during the electrochemical operation. Depending on the quantities accessible via experiments and theory, two different hypotheses are regularly invoked to explain the findings in the range of 0% (graphite) to 100% (LiC$_6$) state of charge (SOC): the staging and the domain model. The lithium intercalation process shows evidence of multiple phase transitions in the voltage vs. SOC diagram. The corresponding system configurations are termed “stages” I, II and so forth. In the simple staging model, these correspond directly to the numbers of empty galleries (spaces between graphene sheets) between the fully occupied ones (see Figure \[fig:staging\]).

![Sketch of Li-intercalated graphite in stage I to III configurations [@Smith2017]. Violet spheres represent lithium ions, black lines correspond to graphene sheets. Bottom right: illustration of the domain model [@Daumas1969]. The structure has the same nominal stoichiometry as the structure in stage II (top right). \[fig:staging\]](./staging_compact.png "fig:"){width="0.9\linewidth"}\

In the domain model, these motifs are not assumed to range over meso-/macroscopic dimensions but to form regions of finite lateral extent. Consequently, it is quite clear that different SOC with the same nominal stoichiometry LiC$_x$ will not be configurationally homogeneous, making Li-intercalated graphite a profoundly non-trivial system to address.

In order to effectively connect to experimental studies, a theoretical framework for simulating large-scale and long-duration non-equilibrium processes in the graphite anode, based on kinetic Monte Carlo (kMC) [@Andersen2019] simulations is required. The first step towards this goal is gaining the ability to quickly and accurately calculate diffusion barriers on the fly, which is the primary motivation of this work. This requires the ability to reproduce reliably and accurately the layer distances (ideally of all possible configurations, but predominantly of the dilute, low-saturation stages) and the forces affecting the lithium-ions, while the strains within the graphene layers are of lesser importance.

Large-scale atomistic simulations typically pursue force field approaches [@Duin2001] for those systems where energetics and kinetics are well described within the upper end of the SOC range. However, those approaches are limited when it comes to the entire range of different SOC, from extremely diluted stages to fully concentrated ones. Recently, a Gaussian Approximation Potential (GAP) was reported to be able to describe amorphous carbon well. [@Deringer2017] However, when the latter was later extended to model lithium intercalation, [@Fujikake2018] it became apparent that the insertion of lithium into those host structures requires a non-trivial description of the electrostatic interaction. Contrary to most approaches, including the one presented in this work, Fujikake *et al.* did not treat the full Li-C system, but attempted to model the energy and force differences arising from lithium intercalation separately, and then added them to the carbon GAP. More specifically, their machine learning process (ML) is based on fitting the energy and force differences between identical carbon host structures, but with and without an intercalated lithium atom. However due to the fact that the lithium intercalation energies are significantly larger in magnitude than the electrostatic lithium-lithium interaction energies, they were not able to recover the latter from the data to a satisfactory degree and had to manually add an extra correction term (fitted to DFT) in order to account for those contributions. To avoid similar shortcomings, we rather base our approach on Density Functional Tight Binding (DFTB) [@Elstner1998], a semi-empirical—and thus computationally much cheaper—approximation to Density Functional Theory (DFT), [@Kohn1965] which has been the most common technique for high-accuracy electrochemical simulations for many decades [@Koskinen2009]. However, since DFTB’s speedup is achieved by pre-calculating atomic interactions to avoid calculating them at runtime, this comes at the cost—or rather, initial investment—of pairwise parametrization. As of now, no Li-Li and Li-C DFTB parameters are available. In the following, we combine for the first time the recently developed Particle Swarm Optimization [@Shi1998] parametrization approach as first proposed by Chou [*et al.*]{} [@Chou2016] with a more flexible ML repulsive potential [@Engelmann2018], to obtain finely-tuned parameters for this system—taking advantage of its physics, albeit perhaps at the expense of some transferability. Let us however stress that the parametrization procedure employed here remains completely general, as the system specificity lies entirely in the choice of the training set(s).

The electronic part {#electro}
===================

![image](./bands_lithium.png){width="0.24\linewidth"} ![image](./bands_graphene.png){width="0.24\linewidth"} ![image](./bands_diamond.png){width="0.24\linewidth"} ![image](./bands_LiC6.png){width="0.24\linewidth"}

In DFTB jargon, the so-called “electronic part” includes the semi-empirical band structure and the Coulombic contributions to the total energy of the system. [@Koskinen2009] These depend parametrically on the diagonal elements $\varepsilon$ of the non-interacting Hamiltonian, the Hubbard-$U$ and a confinement potential which is used to cut off the diffuse tails of the basis orbitals. For the free atom, the first two quantities are tabulated for most elements or can be calculated with DFT. However, using the free atom values is an approximation, and the decision whether it is justified must be made carefully on a case to case basis. The confinement potential, on the other hand, is always treated as a parameter. Quadratic [@Stohr2018] and general power-law functional forms [@Wahiduzzaman2013] are commonly used, as well as the Woods-Saxon potential [@Chou2016] (also employed here) which assures a smoother transition to zero in the orbital tails. Each of these parameters needs to be determined for every chemical species present in the system of interest, typically in a non-linear optimization process. In the PSO, each particle then represents a set of parameters ($\left\{\varepsilon\right\}$, $\left\{U\right\}$, and the confinement constants), with which the DFTB interaction is constructed, so that the parametrization can be improved by minimizing a cost function. The central task is thus the definition of a meaningful cost function. Frequently, one uses the weighted sum of an arbitrary number of contributions $f(\sigma^{DFT}, \sigma^{DFTB})$, each providing a measure of the deviation between DFT and DFTB for some system property $\sigma$. Hereby, as we are optimizing the electronic parameters only, the chosen target properties must not depend on repulsion. For our system, we compare the band structures of metallic lithium, graphene and diamond. Additional details on the definition of the corresponding cost function are provided in the SI. Figure \[fig:bandstructures\] shows our resulting band structures. Overall, we recognize decent agreement for all band structures, while some deviations are expected given the minimal basis in DFTB. For example, the pronounced mismatch in the conduction band at the $H$ point in the lithium band structure as well as the incorrectly direct band gap of diamond can be ascribed to this over-simplification in the DFTB model. For the two carbon systems, we see very good qualitative agreement for most regions of the band structures, but notice a small degree of overall compression towards the Fermi level. Given the systematic nature of this imperfection, we speculate that for further improvement, it would probably be necessary to include $U$ and $\varepsilon$ in our parameter space, which would simultaneously increase the dimensionality of the optimization problem. As an additional validation criterion, we examine the charge population (which DFTB provides by default) for the lithium ions in LiC$_6$. Our parametrization produces a value of $0.853\,e$, in agreement with the value of $0.86\,e$ calculated by Krishnan [@Krishnan2013] with Bader charge analysis [@Bader1990]. Given this excellent agreement and also considering the fact that the repulsion potential is capable of quite effectively correcting small imperfections in the electronic part, we decide not to optimize the latter any further in this work—a decision justified in retrospect by the excellent results we present. However, let us still emphasize the opportunity for improvement here, should it eventually become necessary.

The repulsion potential {#rep}
=======================

It is common practice to assume some analytical form for the repulsive potential and fit the chosen functional parameters as to minimize a set of DFT-DFTB force differences [@Koskinen2009]—a protocol easily implemented also for the PSO approach. However, the main limitation and bias results from the choice of said parametrized functional form. It needs to be sufficiently flexible to cover a large space of systems and bonding situations. This typically yields a high dimensional non-linear optimization problem, which might still be insufficient to capture unexpected subtle, yet extremely relevant physical features. We rather adopt the method recently developed by A. Engelmann [@Engelmann2018], which employs Gaussian Process Regression (GPR) [@Rasmussen2006] to create a flexible functional form “on the fly”, while adapting to the physics captured by the training data set, instead of forcing us to guess it *a priori*. In the SI, we give a short introduction to the method and explain the character and effect of the related hyperparameters, referring the reader to Rasmussen [@Rasmussen2006] for the underlying stochastic theory and to Engelmann [@Engelmann2018] for the application to DFTB repulsive potentials. For the global damping, correlation distance, and data noise hyperparameters, we verified (see SI) that a sizeable subspace of the overall hyperparameter-space is appropriate, and choosing pretty much any combination of values within that subspace will produce very similar, correct results. The same is not necessarily true for the cutoff radii of the C-C and the Li-C repulsion. Since the electronic energy contribution is entirely based on just a sum of non-interacting atomic contributions, the repulsion potential has to account for different chemical environments affecting the same type of atom. In a GPR setting it is therefore of paramount importance to sample a sufficiently large set of training data which covers all interatomic distance ranges and chemical environments relevant for a faithful representation of the system studied. Ideally, it should also be ascertained that the model quality is stable w.r.t. the explicit choice of hyperparameters such as the cutoff radii.

The training data {#training}
-----------------

![Interlayer distances for graphite (grey), LiC$_{12}$ (SOC 50%, grey-purple) and LiC$_6$ (SOC 100%, bright purple) as a function of C-C repulsion cutoff trained. Note that for LiC$_{12}$, there are two different layer distances to consider: one for the empty gallery and one for the full gallery. Here, we plot the average of the two. The dashed lines show the experimental layer distances we aim to reproduce (Sources: Trucano [*et al.*]{} [@Trucano1975] (graphite), Vadlamani [*et al.*]{} [@Vadlamani2014] (LiC$_{12}$ and LiC$_6$). The green coloured area represents the range within which we are satisfied with the performance.[]{data-label="fig:vsCC"}](./all_vs_CC_colors.png){width="0.9\linewidth"}

In terms of DFT functional, our starting point is PBE [@Perdew1996], which has been used by the majority of researchers working on intercalation phenomena and is known to describe LiC$_6$ well. However, it does not reproduce the dispersive interaction between graphene sheets. In order to address this, we finally (see “Set 3” below) combine the reference PBE calculation with a Many Body Dispersion (MBD) treatment and the DFTB model with a computationally cheap Lennard Jones (LJ) [@Zhechkov2005] dispersion correction [@Rappe1992]. The rationale for this choice is that PBE should reproduce galleries containing many lithium atoms correctly and LJ-dispersion should predict empty galleries well, while not interfering too much with the PBE-description of the concentrated ones. However, it is unclear, how this interaction shapes out for intermediate, dilute lithium stoichiometries. During our investigations, we find that this approach works somewhat decently, but needs some controlled adjustments (vide infra) in order to produce truly satisfactory results.

As a first guess, we construct a set of training structures (Set 1) which consists of a balanced mix of Li$_n$C$_{36}$ super-cells ($n \in (0, 1, ..., 6)$), in order to represent the entire range of charging states (exemplary structures are shown in the SI). Additionally, those structures are rattled (each atom randomly displaced), as well as compressed or expanded. This procedure yields a smooth distribution of bond lengths and forces. We then train a GPR repulsion potential by matching DFTB against PBE forces for this structural ensemble, aiming at a first, mostly transferable model. The standard LJ DFTB correction is subsequently applied on top of this parametrized DFTB model. With this approach, we are able to find parametrizations that reproduce all layer distances (of graphite, LiC$_{12}$ and of LiC$_6$) correctly, albeit not for a stable range of all parameters (in particular the Li-C cutoff, see below). As shown in Figure \[fig:vsCC\], the choice of cutoff radius for the C-C repulsion potential does not have a major influence on the layer-distances for quite a large range of values. In fact, the point at which the predictions stop being accurate can be identified as approximately the experimental values for the interlayer distances. Going beyond that with the cutoff radius essentially corresponds to including interlayer interactions in the potential fit, mixing their description with the intralayer covalent bonds. Thus, the restriction of the cutoff radius we find here is physically motivated by the range separation of the interactions that characterize our system: as the 2nd next neighbour distance in a relaxed graphene sheet is around $2.45$ Å and the layer distance is $3.35$ Å, the cutoff range defined by the (smallest) plateau in Figure \[fig:vsCC\] represents a sweet spot where the GPR learns 2nd next neighbour interactions but does not yet (mistakenly) take any interlayer interactions (even in the compressed structures) into account in the repulsion potential. In light of these findings, we select the cutoff value $2.6$ Å for the C-C-repulsion potential. Indeed, we did not encounter any reason to change this selection during the entirety of this work (despite rigorously testing it for each of the training data sets).

However, with this first training set we do not obtain an equally stable plateau as a function of the Li-C repulsive cutoff (see SI). Furthermore, we discover that the quite strongly distorted graphite planes in these structures lead to large forces compared with those acting on the intercalated lithium-ions hindering the performance in lithium-forces prediction. We tackle the second problem first: while the rattled, scaled structures in Set 1 cover a sufficiently large range of bond lengths, they only account for configurations with the lithium-ions sitting over the centre of a graphite ring, *i.e.* in a local energy minimum. We recognize this as the reason for the comparably small lithium-forces. In order to balance out this structural bias, we calculate a number of transition paths for lithium diffusion processes using a Nudged Elastic Band (NEB) method [@Henkelman2000; @Henkelman2000a]. Exemplary structures can be found in the SI. Now, we are able to extract structures from these trajectories, in which the lithium ions are subject to stronger forces commensurable with the graphite-layers. For our second training set (Set 2), we replace around $45\%$ of the rattled and scaled structures with those transition path geometries. By this measure, we are able to improve the accuracy for predicting forces on Li-ions significantly, without sacrificing the description of the graphite layers. However, while we do observe a plateau for the resulting layer distances with respect to the Li-C cutoff, the interlayer distances are not reproduced equally well as in Figure \[fig:vsCC\] for Set 1 (see Figure \[fig:vsLiC\], yellow area), while only the LiC$_{12}$ interlayer distances assume correct values (see Figure \[fig:vsLiC\], green areas), yet outside the plateau.

![Interlayer distances for LiC$_{12}$ (SOC 50%, grey-purple) and LiC$_{6}$ (SOC 100%, bright purple) as a function of Li-C repulsion cutoff, with a fixed C-C cutoff set to 2.6 Å. The repulsion was trained on a set analogous to Set 1 ([*cf.*]{} text), where  45% of the structures were replaced by geometries randomly extracted from intra-layer Li diffusion paths. For LiC$_{12}$, the plotted interlayer distance is the average between the values for the filled and the empty gallery. The dashed lines show the experimental layer distances. The yellow coloured area represents the range within which the results are stable, however at a wrong value.[]{data-label="fig:vsLiC"}](./all_vs_LiC_set2.png){width="0.9\linewidth"}

![Interlayer distances for LiC$_{12}$ (SOC 50%, grey-purple) and LiC$_{6}$ (SOC 100%, bright purple) as a function of Li-C repulsion cutoff, with a fixed C-C cutoff set to 2.6 Å. The repulsion was trained on a set analogous to Set 2 ([*cf.*]{} text), where  55% of the structures were replaced by geometries with MBD-corrected forces. For LiC$_{12}$, the plotted interlayer distance is the average between the values for the filled and the empty gallery. The dashed lines show the experimental layer distances. []{data-label="fig:vsLiC_set3"}](./all_vs_LiC_set3.png){width="0.9\linewidth"}

![image](./landscapes.png){width="0.9\linewidth"}

This behaviour suggests that our problem here does not lie in the choice of the training set, but rather in the treatment of long-ranged interactions.

Let us consider the underlying predicament: so far, the DFTB-part of the force residues used for the ML process is calculated without LJ dispersion correction. We then construct the repulsion potential with the purpose of making those DFTB calculations match references based on PBE-DFT, which reliably predicts layer distances for LiC$_6$. By then using LJ (required to obtain the correct empty layer distance in graphite) in our actual DFTB calculations (after the parametrization process), we cause the aforementioned offset for highly lithiated compounds. Using LJ already for the force-residue calculations during the ML seems like the obvious solution to this problem. However, this presents a new issue in the lower-saturation range (Li$_x$C$_6$, $x<0.5$). There, we previously fitted the repulsion to PBE-DFT references, which are not correct in that range without dispersion correction. The resulting DFTB forces are then shifted by LJ towards the correct value (as is indicated by the quite decent results for LiC$_{12}$ with Set 2). But after the modification, we would then fit the *final* DFTB forces (that result after applying the LJ) to the (incorrect) PBE-DFT references, thus improving our performance for highly saturated system states, but ruining it for dilute ones, by effectively double counting dispersive contributions. It becomes apparent that in order to make this approach work, we need to utilize dispersion corrected DFT reference forces which are also correct for low saturation states and, at the same time, compatible with the computationally cheap DFTB-LJ correction.

Our ansatz is that we can essentially—to a degree—encode the difference between the LJ dispersion and the “true” dispersion into the repulsion potential. At this point we stress that *ideally*, both the true, non-local exchange correlation functional in DFT and an ideal repulsion energy in DFTB would already encompass all dispersion effects, and it is solely due to approximations in the derivations, [*e.g.*]{} of GGAs, that they do not in these models. Therefore, rather than mixing our repulsion potential with something fundamentally different (which would be physically questionable), what we do here simply corresponds to partially adding a contribution back in, that should have been there in the first place. To our knowledge, the currently best way to calculate dispersion corrected lithium intercalated graphite, with correct layer distances predicted for the entire saturation range, is the MBD correction [@Tkatchenko2012]. This method is computationally rather expensive, but since we only need to run DFT calculations for our training data set, which is very limited in size, this is not vital to us. We do realize that this approach most likely comes with some cost in terms of transferability. In order to retain as much of it as possible, we choose not to replace *all* force residues, but only $\approx50\%$, which proves sufficient to demonstrate the effectiveness of the presented method in a general way. Nonetheless, further investigating the effect this percentage has on the performance is certainly a task that should be tackled in the future. Of course, alternatively to our approach, it is possible to simply apply the MBD correction scheme directly to our DFTB calculations. However, doing so would cost us one to two orders of magnitude in speed, as MBD then becomes the dominating step in terms of computation time. Exactly as we had hoped for, we have succeeded at shifting the predicted interlayer distances (within the stable Li-C cutoff plateau) into the very close proximity of the experimental reference values for both LiC$_6$ and LiC$_{12}$ (Figure \[fig:vsLiC\_set3\]). Especially the excellent results for the stage II compound LiC$_{12}$ show that our parametrization is now able to handle *both* mainly ionic concentrated *and* mainly dispersive dilute layers to a satisfactory degree. In Figure \[2Dlandscape\], we illustrate the effect our modification has on the repulsion potential landscape for a wide range of Li-C cutoff radii. First (and most notably), we have moved and solidified the local minimum related to the next-neighbour lithium-carbon interaction (see bottom right). For the Set 2 and Set 3 potentials, the minima (blue and black dashed lines) are located at atomic distances of $2.41$ Å and $2.35$ Å respectively, which correspond to LiC$_6$ interlayer distances of $3.83$ Å and $3.67$ Å, the exact values which *do*, in fact, result from the relaxation of those structures, using the two repulsion potentials respectively. The 2D maps (top) show that this behaviour is apparent for an entire range of cutoff radii, thus ruling out the possibility that the fit is only accidentally correct (as it happens, [*e.g.*]{}, for Set 2, see Figure \[fig:vsLiC\]). In bottom left, we can also clearly see the upper ($\sim 5.8$Å) and lower ($\sim 4.5$Å) boundaries for the cutoff radius, beyond which the physicality of the model falls apart. They define exactly the range within which we find the stable cutoff dependency plateau, which is now at the correct numerical value as shown in Figure \[fig:vsLiC\_set3\]. We may identify the upper boundary at $5.8$ Å (as further discussed in the SI), as the distance between a lithium ion and the second closest graphene sheet, which is an intuitively plausible limitation. It is less obvious, though, to assign a clear physical meaning to the lower bound at $4.5$ Å, as it cannot be directly related to any particular structural feature of LiC$_6$. The most likely cause, we believe, is that the cosine-shaped cutoff function employed in the GPR framework starts cutting off physically relevant details from the repulsion potential below that. A physically motivated lower bound may be identified by evaluating the mean absolute forces acting on lithium as a function of Li-C cutoff, shown in Figure \[Fig:maf\].

![Mean absolute forces acting on Li for validation set??? compared to DFT reference calculations (black dashed line)[]{data-label="Fig:maf"}](./rmsd_forces_set3.png){width="0.9\linewidth"}

Overall, we now observe two separate Li-C cutoff plateaus: between approximately 4.3 Å and 5.7 Å, we obtain accurate layer distances (Figure \[fig:vsLiC\_set3\]), while for radii above roughly 5.0 Å, our predictions for forces and transition energies are correct (Figure \[Fig:maf\]). This duality can very simply be explained by the fact that the first property is mostly a z-direction phenomenon (and interactions with the second closest graphene sheet limit the physicality of our model), while the other takes place almost exclusively in the x-y-plane, where no such limitation applies. Given this difference in fundamental nature, it is very plausible to trust both these plateaus. Thus, their overlap (5.0–5.7 Å) defines the region within which any value of the Li-C cutoff radius produces an almost identical parametrization that performs very well, for all our benchmark criteria, in a stable and trustworthy manner.

Results: interlayer distances and diffusion barriers
====================================================

[l l | r r r r r r ]{}\

&compound &experiment &DFTB &DFT in [@Krishnan2013] &filled gallery &empty gallery &barrier\
C$_6$ &– &$3.355$ Å &$+46$ mÅ &$+62$ mÅ &– & $+46$ mÅ\
LiC$_6$ &stage I &$3.687$ Å &$-12$ mÅ &$+56$ mÅ &$-12$ mÅ &– &$351$ meV\
LiC$_{12}$ &stage II &$3.511$ Å &$+16$ mÅ &$-16$ mÅ &$+112$\* mÅ &$-79$\* mÅ &$401$ meV\
LiC$_{18}$ &stage III &$3.470$ Å &$-30$ mÅ&$+173$ mÅ &$+196$\* mÅ &$-50$\* mÅ &$393$ meV\
\

\

Table \[table1\] reports some resulting inter-layer distances and diffusion barriers based on our new DFTB parametrization in Table \[table1\], compared with experimentally determined values, as well as previous theoretical findings.

Furthermore, we draw qualitative conclusions from these results and summarize their implications on the intercalation mechanism. As a quick reminder, stages I, II and III correspond to every, every other and every third gallery being filled (to any degree) with lithium. Additionally, we describe the concentration of the intercalant in a filled gallery as dilute (low) or concentrated (high), thus allowing for a simple classification of fundamentally different compounds.

Here, we take only concentrated stages into consideration. For all calculations, we chose an Li-C cutoff radius of $5.5$ Å, following the findings discussed above. As Table \[table1\] clearly illustrates, we systematically outperform the method by Krishnan *et al.* [@Krishnan2013]—in terms of accuracy—for every structure they provide comparison for. This is especially remarkable considering the fact that they used full GGA-DFT with dispersion corrections in post-processing, which is the current state-of-the-art approach, as well as significantly more computationally expensive than our method.\
Subsequently, we investigate intra-layer next-neighbour diffusion barriers and compare our results to recent experimental findings from [@Umegaki2017] (based on muon spin relaxation spectroscopy). Our calculations yield purely *microscopic* results within 50 meV from each other for all three relevant compounds, as is shown in Table \[table1\]. The deviations between them can be attributed to slight differences in the filled-layer spacing of the different structures.\
In contrast, the experimentally determined *active* barriers of 270 meV for LiC$_6$ and 170 meV for LiC$_{12}$ show a strong dependency on the systems stage. We believe this difference to be caused by correlation effects. Capturing those using kinetic Monte Carlo simulation is something we intend to do in the near future.

Conclusions and outlook
=======================

In this work, we put forward—for the first time combining particle swarm (i.e. PSO) and machine learning [@Engelmann2018] (i.e. GAP) approaches for this task—a well-performing DFTB-parametrization for lithium intercalated graphite which is capable of very accurately reproducing various structural properties and qualitative trends relating to the intercalation mechanism for a wide variety of Li$_x$C$_{36}$ compounds. In the course of this process, we believe to have shown that Density Functional Tight Binding (DFTB) is a superior approach for modelling intercalation compared with force field methods (*e.g.* the GAP by [@Fujikake2018] requires a manual correction term for lithium-lithium interactions which our method does not). Furthermore, we share key details and choices along this process and thus provide guidance for similar endeavours in the future.
