Accessing Structural, Electronic, Transport and Mesoscale Properties of Li-GICs via a Complete DFTB Model with Machine-Learned Repulsion Potential

Lithium-graphite intercalation compounds (Li-GICs) are the most popular anode material for modern lithium-ion batteries and have been subject to numerous studies—both experimental and theoretical. However, the system is still far from being consistently understood in detail across the full range of state of charge (SOC). The performance of approaches based on density functional theory (DFT) varies greatly depending on the choice of functional, and their computational cost is far too high for the large supercells necessary to study dilute and non-equilibrium configurations which are of paramount importance for understanding a complete charging cycle. On the other hand, cheap machine learning methods have made some progress in predicting, e.g., formation energetics, but fail to provide the full picture, including electrostatics and migration barriers. Following up on our previous work, we deliver on the promise of providing a complete and affordable simulation framework for Li-GICs. It is based on density functional tight binding (DFTB), which is fitted to dispersion-corrected DFT data using Gaussian process regression (GPR). In this work, we added the previously neglected lithium–lithium repulsion potential and extend the training set to include superdense Li-GICs (LiC6−x; x>0) and lithium metal, allowing for the investigation of dendrite formation, next-generation modified GIC anodes, and non-equilibrium states during fast charging processes in the future. For an extended range of structural and energetic properties—layer spacing, bond lengths, formation energies and migration barriers—our method compares favorably with experimental results and with state-of-the-art dispersion-corrected DFT at a fraction of the computational cost. We make use of this by investigating some larger-scale system properties—long range Li–Li interactions, dielectric constants and domain-formation—proving our method’s capability to bring to light new insights into the Li-GIC system and bridge the gap between DFT and meso-scale methods such as cluster expansions and kinetic Monte Carlo simulations.


Introduction
Lithium-graphite intercalation compounds (Li-GICs) are the primary anode material for commercial Li-ion batteries with a market share of 98% [1] due to their good volumetric and gravimetric capacities, long cycle life, abundant availability, and low cost. Despite investigations into alternatives such as lithium-metal anodes, graphite and modified-graphite compounds will not be replaced in the foreseeable future, as important EV manufacturers, material suppliers and cell producers have recently announced that graphite-containing composites will mark the state of the art for next-generation lithium-ion batteries [2].
Additional details on how the potential shapes changed between [26] and the present work are provided in the Supplementary Information. DFT calculations, serving as a reference for the DFTB fit, were performed with the all-electron framework FHI-aims [29] with light settings and default tier-2 basis sets, using the PBE exchange-correlation functional [30]. For dispersion correction, the MBD approach was chosen [31,32]. MD simulations for generating the training set and validation structures were performed in the NVE ensemble at 300 K and 1000 K using the LAMMPS code [33] with the embedded atom method (EAM) potential for alkali metals developed by Nichol and Ackland [34].
Geometries were constructed and analyzed by means of the atomic simulation environment (ASE [35]) which we also used as a base framework for all force and energy calculations, structure relaxations (specifically using the BFGS algorithm as an optimizer [36]), and barrier calculations. For the latter, we employed the nudged elastic band (NEB) [37,38] algorithm with the FIRE-optimizer [39] and climbing image switched on.
For all DFTB calculations, we used a well-converged k-point density of at least 0.1/Å. The SCC-tolerance is 10 −7 . We employed Fermi filling with a Fermi temperature of 300 K, as well as a Broyden mixer [40] for convergence acceleration with a mixing parameter of 0.5.
All of these settings were tested with regard to convergence for the whole range of SOC. As described in [26], our parametrization is meant to be used with the Leonnard-Jones dispersion correction [41] switched on.

Experimental
Open circuit voltage (OCV) curves were recorded to compare the simulation with real measured values. For this purpose, cells were built with graphite against lithium as well as highly oriented pyrolytic graphite (HOPG) against lithium. The cells from EL-Cell (ECC-Std), comparable to button cells, were used as the housing. A Whatman GF\D was used as the separator and EC:DMC 1:1 with 1 mol/L LiPF6 from Sigma-Aldrich was used as the electrolyte. The graphite electrodes were coated on copper current collectors, whereas the HOPG was used without a current collector. The electron conductivity was sufficient due to its low current rate. The graphite was used in 18 mm blanks, whereas the HOPG was cut into narrow strips with a width of approximately 2 mm to ensure the highest possible surface-to-volume ratio. The ions can only intercalate into the HOPG from the cut edges and not through the surface. The graphite cells were initialized with a current rate of C/10 and the HOPG cells were cyclized with a current rate of C/30.
All considered benchmark methods performed well in reproducing the in-plane lattice parameter a, which is governed by covalent C-C bonds. However, even state-of-the-art dispersion-corrected DFT functionals (except for [14]) struggle with predicting the out-ofplane lattice parameter c. This is due to the fact that the latter is governed by van der Waals interactions, which are still notoriously difficult to account for despite considerable effort in creating various correction schemes [31,[44][45][46] for DFT. The recent machine learning approach [10] overestimates c by an even larger margin.
GPrep-DFTB results are very close to the experimental references-closer than even the majority of DFT approaches-proving that the method is very capable of treating both covalent C-C bonds and van der Waals interactions in graphite.
When lithium intercalates into graphite, the graphene sheets shift from AB-stacking to AA-stacking somewhere between 5% and 15% SOC, and the interlayer distance expands from~3.36 Å to~3.62-3.7 Å (depending on the SOC of adjacent galleries) for the full gallery. Empty galleries adjacent to filled galleries also slightly expand, due to the extra charge transferred to the graphene sheet from the intercalated Li-ions, making the overall increase in the average z direction lattice parameter non-linear.
It is generally accepted that Li-ions do not evenly distribute throughout the entire GIC, but tend to arrange themselves in fully filled domains and empty domains [7,47] (see Figure 1), leading to a local staging behavior with the staging number (n = 1, 2, . . .) indicating that n − 1 galleries are empty between each pair of filled galleries only within that limited region. For our second benchmark, we calculated the average interlayer distances depending on the stage n = 1, . . . , 9 of stoichiometry LiC 6n (corresponding to LiC 6 , LiC 12 , LiC 18 , LiC 24 and higher) with the GPrep-DFTB and compared with experimental and DFT references, where available ( Figure 2 and Table 2).
In order to be able to directly compare with DFT, this set of unit cells was constructed with global staging (Figure 1, left) and not according to the domain model, which would render them far too big for DFT. Because of that, it is not obvious whether AAor ABstacking should be assumed for the empty layers. In real samples, which are large of scale and governed by the domain model ( Figure 1, right), it is probable that empty parts of galleries also exhibit AA-stacking, because they are forced into that configuration by adjacent filled domains within that same gallery and because they are not truly empty, either. However, in an idealized system, without factoring this in (Figure 1, left), ABstacking of the empty galleries (as in pure graphite) is also conceivable. Therefore (and because it is unclear which stacking order has been assumed in the DFT reference [13]), we provide predictions for both as well as a prediction area (light blue). Filled galleries are always in AA-stacking. For the empty galleries, we predict that interlayer distances remain mostly constant throughout all stages with a slight increase in stage 2, due to the additional electrostatic repulsion caused by the charge transfer from Li-ions in the adjacent filled galleries to the carbon sheets. For the filled galleries, a constant interlayer distance can be observed for stages 4 and higher, whereas stages 3, 2, and 1 progressively show increased interlayer distances, which we attribute again to the electrostatic repulsion of the increasing charge density. Based on these findings, a simple building block model that assumes invariant interlayer distances proves sufficient to describe the system's behavior in the high-stage (n > 3) limit, where filled galleries are too far apart to interact in any way ( Figure 2, left: 'fit AA' and 'fit AB'). Only for stages 1 and 2 can the increased filled interlayer distance cause a significant deviation from this model.
Experimental results [5,9] agree very well with our GPrep-DFTB predictions for stages 1, 2, and 3. The DFT-based Ising model by [13] is also accurate for stages 1 and 2, but maintains an overly steep slope for higher stages 3 and 4, which-if continued-would clearly converge towards wrong asymptotic behavior. AA and AB signifies AA-stacking and AB-stacking assumed for the empty graphite layers. Where available, experimental and theoretical references are also shown. The curves 'fit AA' and 'fit AB' correspond to a simple building block model with just two fixed interlayer distances for empty and filled galleries, respectively. For stages 3 and higher, our predictions adhere closely to that model. For stages 1 and 2, the interlayer spacing of the filled galleries was expanded due to the additional charge transfer. (Right): Interlayer spacing of the full galleries, as well as empty galleries in AA and AB stacking, as a function of the stage (calculated with GPrep-DFTB).

Diffusion Barriers
Diffusion barriers for ion transport are among the most interesting properties of mixed ion-electron conductor (MIEC) materials such as Li-GICs. They are the crucial input parameter of any kinetic Monte Carlo simulation [48] that aims to predict large-scale phenomena such as phase transitions between stages or non-equilibrium configurations during fast charging. They are also quite difficult to reliably calculate, since they are closely linked to the interlayer distance, which, as previously pointed out, is still hard to predict, even with state-of-the-art DFT, especially for dilute, low SOC configurations.
In order to rigorously investigate the transport properties of Li-GICs, we constructed a variety of structures based on 2-and 3-layered supercells with 36 and 48 carbon atoms, respectively. This allows us to consider different staging and Li-stacking orders for equivalent stoichiometries. We fully relaxed all structures, extracted the interlayer distances, and calculated the energy barriers (Table 3) for exemplary bridge-path diffusion processes (connecting the next-neighbor Li positions) using the NEB method (see Supplementary  Information). for the exact initial and final states of the bridge path NEBs, as well as the predicted barrier diagrams. Analyzing the interlayer distances first, there are multiple trends to observe. First of all, structures with an AαAα order (Greek letters indicating the Li layer order) consistently relax to a smaller interlayer distance than structures with an AαAβ stacking, implying that the former may be more favorable. In terms of the total energy per unit of 6 carbon atoms, however, we see virtually no difference (AαAβ: −297.279 eV/6C, AαAα: −297.253 eV/6C) for both structures fully relaxed in terms of cell parameters and all atom positions. This deviation of E(AαAα) − E(AαAβ) = 26 meV/6C is on the order of k B T at room temperature, implying that at ambient conditions, no clear distinction between the Li orderings can be made and experiments would probably see a mixture of both. For reference, DFT (PBE + D3), which we consider trustworthy at least for high SOC compounds, predicts a deviation of E(AαAα) − E(AαAβ) = −14 meV/6C. It is necessary to recognize the difference in sign here, but since both values are within k B T at room temperature and at the limit of DFT errors, we do not believe that this constitutes a relevant difference.
Furthermore, for stoichiometries which can either be arranged as dilute stage 1 or dilute stage 2 (Li 2 C 36 ) geometries, stage 2 has the lower interlayer distance and is therefore favored, which is in line with the agreed upon theory of staging and domain formation [47]. This is due to the fact that the z direction expansion of a gallery does not vary linearly with the filling factor. The total expansion from the 0% filling to 100% filling is 0.315 Å for AαAα-stacked stage 1 compounds. Filling empty layers by 33% (stage 1-Li 2 C 36 ) already expands the interlayer distance by 0.165 Å, which is 52% of the total expansion. At 66% filling (stage 1-Li 4 C 36 ), the interlayer distance is expanded by 0.252 Å, which is 80% of the total expansion.
In terms of the diffusion barriers, experimental sources vary quite a lot. This is due to the fact that the measuring technique, as well as additional factors such as temperature, play a role. Furthermore, if total diffusivities are measured, it is difficult to separate those into energetic and kinetic contributions. Langer et al. [21] measured a barrier of 0.55 eV (LiC 6 ) by means of Lithium nuclear magnetic resonance (Li-NMR), Freilander et al. [22] report 0.6 eV (LiC 6 ) and 1.0 eV (LiC 12 ) using beta-NMR and Magerl et al. [23] find 1.0 eV (LiC 6 ) by means of the quasielectric neutron scattering (QENS) at T > 630 K.
On the theoretical (DFT) side of things, the revPBE-D3-BJ approach by Thinius et al. [14] (which has proven very accurate for structural parameters and formation energies) predicts the barriers of 0.42 eV (LiC 6 ) and 0.47 eV (LiC 12 ). Toyoura et al. [24] reported 0.3 eV (LiC 6 ) and 0.49 eV (empty gallery) by means of DFT (LDA) and Persson et al. [25] predict 0.283 eV (LiC 6 ) and 0.297 eV (LiC 12 ), using DFT (GGA) with the interlayer distances fixed at experimental values. Even though there is some variation in the absolute numbers, both theoretical and experimental studies are in general agreement on the ordering of the barrier heights: LiC 6 < LiC 12 < empty gallery. This corresponds to an inverse dependency on the interlayer distance of the respective gallery.
Using GPrep-DFTB, we predict barriers of 0.404 eV (LiC 6 ), 0.426 eV (LiC 12 ) and 0.504 eV (empty gallery). These results capture the same previously explained qualitative trend as the references. This also holds true for other configurations that had not been investigated before, such as AαAβ-stacked and stage 3 structures. Quantitatively, our results are in particularly good agreement with the revPBE-D3-BJ approach [14].

Formation Energetics
The intercalation energies of Li-ions entering the GIC at different states of charge are a crucial measure for predicting the most stable configurations throughout the charging process and the phase transitions between those. Therefore, as a third benchmark, we calculated the formation energies of LiC 6 , LiC 12 and LiC 18 , and compared them with various experimental and theoretical studies ( Table 4).
We calculate the intercalation energies per lithium atom (or, equivalently, per formula unit), as where n is the stage number and E(·) is the DFTB total energy. These are directly comparable to the corresponding literature values obtained by DFT calculations. The latter thus do not include any finite temperature effects. On the other hand, experimental values are formation-free energies or enthalpies. Additionally, the calorimetric reference [20] is taken at elevated temperature (455 K) and with liquid lithium as precursor, rather than at room temperature with a solid lithium electrode. As a final note, the calculated values correspond to infinite phases of perfect stage n stoichiometry, whereas the true compounds at the corresponding stoichiometries are a mixture of domains of yet unknown structural details. Consequently, the values given within the scope of this parametrization study are not yet intended to be quantitatively comparable to the experiment, but one may still identify qualitative trends, just like with regular DFT. The DFTB model opens up the way to forthcoming more realistic, quantitative simulations.  [20]; DFT references: (d) revPBE-D3-BJ [14]; (e) optB88-vdW [12]; ( f ) revPBE-vdW [12]; (g) GGA-PP [11]; * These values are intercalation free energies ∆G int . The value for LiC 18 is not directly given in the paper but can be consistently extracted using the same formula the authors used for LiC 12 .
For the formation energies, the experimental studies [6,19,20] do not agree as closely, as for the structural parameters, but they do at least provide the same ordering for LiC 6 and LiC 12 and LiC 18 . Interestingly, the values extracted from open circuit voltage (OCV) measurements only agree in the ordering if the energies are taken per formula unit. Normalizing the energies per graphite unit, LiC 12 would have a less negative formation energy than LiC 6 . As the OCV curves in [6,19] agree with each other, we attribute the mismatch to different methods for extracting the formation energies from the OCV curve. We also note here that our structural models closely correspond to highly oriented pyrolytic graphite (HOPG), while all the referenced OCV curves were taken with different forms of graphitic carbon. In order to make sure that the deviation was negligible, we measured the OCV curve for HOPG. Due to the small insertion surface for Li-ions, the characteristic voltage plateaus are not as visible. However, in the regions corresponding to the phase transitions between LiC 18 to LiC 12 and LiC 12 to LiC 6 , the measured HOPG curve matches within 0.01 V with the references. Our measured OCV curve is shown in the Supplementary Information for both HOPG and the coin cell.
Similarly, as for the C-C interlayer distance, the different DFT functionals vary significantly in their performance predicting the formation energetics of Lithium-GICs. Compared with the experimental studies, revPBE-D3-BJ [14] proves to be best, just as it did previously for the C-C interlayer distance. According to [12], revPBE-vdW correctly predicts the phase transitions between stage 1 and stage 2 qualitatively (even though it strongly underestimates the formation energies), whereas optB88-vdW does not. For revPBE-D3-BJ, we do not have this information.
Overall, the formation energy for LiC 6 tends to be consistent across experimental measurements and most computed references. Our results with GPrep-DFTB are equally accurate for LiC 6 , while for both LiC 12 and LiC 18 , we obtain more negative formation energies than the majority of references (with the exclusion of [6]). However, this is not necessarily a pitfall, considering that the finite temperature contributions are not included. The effect of the latter is generally nontrivial; in particular, the entropy variation in Li-GICs is negative for the largest part of the state of charge range. Given the overall uncertainty in both experimental and computed references, we leave this question open for further investigation and adjustments to the parametrization, if needed. We note in passing that potential refinements to the GPrep-DFTB parametrization are possible with little effort, by simply retraining the repulsion potential with additional training data and/or finely tuned hyperparameters. As a perspective, we intend to use this parametrization to train a cluster expansion similarly to [10,49], in order to perform free energy sampling and calculate the OCV curve. If that agrees with the measurements, then the non-perfect energetics of single ideal geometries is only a minor setback.

Long-Range Interactions
Having successfully benchmarked our GPrep-DFTB approach against a variety of comparatively small-scale properties, we proceeded towards calculating some larger-scale properties which are out of reach of DFT (at least at a reasonable computational cost). First, we want to investigate the long-range in-plane interaction between two Li-ions within Li-GIC. In order to do that, we constructed a supercell with 218 carbon atoms and two layers to eliminate any periodic next-layer interactions and allow for the bulging of the graphene sheets. We then exhaustively performed 47 full structure relaxations of all symmetry-inequivalent local minima and maxima for two Li atoms within a single layer in that supercell, as well as 41 five-image NEBs for the diffusion processes between each adjacent pair of local minima. With our method, all of this is possible within days and on a regular workstation. This leaves us with 170 data-points, which we use to fit a 2D potential energy landscape for the whole supercell ( Figure 3a) and also to plot the Li-Li interactions as a function of the Li-Li distances (Figure 3b,c).
As our results clearly show, the in-plane Li-Li interactions are governed by Coulombic repulsion. Even quantitatively, our predictions agree very well with the approach of Pande et al. [13] (BEEF-vdW-DFT + Ising model). However, they were only able to provide four data points which are local minima and therefore quite cheap to calculate, whereas we can also predict transition states (which require NEB calculations and are much more computationally expensive because of that). This proves that GPrep-DFTB is capable of very accurately capturing the Li-GIC system's electrostatic properties, and of doing so for vastly more and larger supercells than DFT.  Based on this, it is then possible to extract the slope from Figure 3c, which, via the relation: gives us access to the relative dielectric constant ε r of the system. We note, however, that this is the effective dielectric response experienced by Li-ions within the Li-GIC and not the macroscopic dielectric constant of the GIC including the contribution due to the Li ion motion. By means of linear regression, we obtain a slope of 0.996 ± 0.015 eVÅ. This leads to ε r /Z 2 = 14.46 ± 0.22. For an assumed partial charge Z of 0.8 to 0.9 for the Li-ions, which is in line with the charge analysis performed by Rana et al. [50], we then predict a relative dielectric constant of 9.1-11.9. Expanding on this in future work, we will be able to, for the first time, determine the dielectric constant of Li-GIC as a function of the state of charge, which is an important input parameter for kinetic Monte Carlo simulations of charge carriers.

Domains vs. Dilute
While the general truth of the Daumas-Heróld domain model [7] has been widely accepted and supported by both theoretical [51,52] and experimental [53,54] studies, qualitative details such as domain sizes and shapes, dependencies on the charging speed and other dynamic factors have not been understood to a sufficient degree. As pointed out in [55], the formation of domains is at least partially responsible for the wide range of Li-ion diffusivities reported from the experiment (10 −6 − 10 −14 cm 2 s [25]), and therefore of crucial importance for understanding the overall behavior of Li-GICs, especially when exposed to the rapidly growing charging speeds that are necessary today.
In order to provide a further demonstration of our method's potential to bring forward a new understanding of these phenomena in the future, we constructed four supercells with roughly 600 atoms each-two of them in a dilutely spaced configuration (with the Li-ions spread evenly across the filled layers, right column in Figure 4) and two of them with a configuration according to the Daumas-Heróld domain model (Figure 4, left column)-and performed full structure relaxations on each of them. Note that what we show here is one possible realization of the dilute LiC 12 and LiC 24 , but there are necessarily many other disordered realizations with very similar total energies and any experiment would likely see a mix of these. Our results show that both in the stage 1 (upper row) and the stage 2 (bottom row) compound, the domain-structure expanded to a significantly smaller overall interlayer distance, meaning it is favored compared to the evenly spaced one. This agrees with the Daumas-Heróld domain model [7]. Furthermore, we were also capable of extracting local interlayer distances for different areas of the structure. As shown in Figure 4, the difference in interlayer distance between the empty and filled domains is larger than 0.2 Å, which is much larger than any residual differences in the computed values above. According to both DFT [56] and our own results (Table 3), this difference corresponds to a difference in diffusion barriers of approximately 0.1 eV or 25%. Given this direct dependency, GPrep-DFTB could, for example, be used to reliably predict local diffusivities in large structures without the need to perform costly NEB calculations. Additionally, we believe that the capacity of GPrep-DFTB for fully relaxing large unit-cells with a multitude of Li-ion distributions has great potential for building more diverse and well-rounded training-sets for, e.g., lattice gas expansions or machine-learning models, than DFT could, also including the 'empty' or 'almost empty' regions of the Daumas-Heróld domain model, for which we provided an extensive model in Figure 3. This makes GPrep-DFTB a crucial new bridge between the atomistic scale and the macroscopic scale of Li-GIC modeling.

Conclusions
The structural, energetic and electronic properties of Li-GICs were theoretically investigated with our DFTB parametrization (based on GPR repulsion fitting) and benchmarked against dispersion-corrected DFT calculations and experiments. The calculated lattice parameters of graphite (a = 2.476 Å, c = 6.746 Å) agree better with experiments than most DFT approaches. For stages 1 through 4, our method correctly predicts the nonlinear nature of the increase in the interlayer distance in LiC x upon intercalation, not only qualitatively but quantitively as well. LiC 6 relaxes to an interlayer distance of 3.682 Å. The calculated formation energies of −0.14 eV (LiC 6 ), −0.40 eV (LiC 12 ) and −0.42 eV (LiC 18 ) per formula unit slightly overestimate the experimental results, but are within the range of DFT predictions. We expect future calculations which include entropy effects to be even more accurate. The calculated diffusion barriers (0.396 eV-0.504 eV, depending on the configuration) show trends supported by accepted theory, such as the Daumas-Heróld domain model, and agree with state-of-the-art DFT studies and experiments. In terms of long-range Li-Li interactions, our model captures the Coulombic nature also discovered by DFT, but is at the same time capable of accessing much larger supercells. Based on these calculations, we predict a dielectric constant for LiC 108 in the range of 9.1-11.9 and recognize the potential of GPrep-DFTB to, for the first time, calculate the dielectric constant of the Li-GIC as a function of the SOC in the near future. Finally, the GPrep-DFTB relaxation of large structures in both dilute and domain-like configurations predicts less expansion of the interlayer distance for the domain structure-again agreeing with previous studies and illustrating the potential of our method for further investigation into the complex and large-scale physics taking place in Li-GICs and for being a new kind of bridge between the atomistic scale and the macroscopic scale of future battery materials.