{'fraction_non_alphanumeric': 0.02317462675480948, 'fraction_numerical': 0.012812894600014856, 'mean_word_length': 4.713346064078082, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Robust high-fidelity DFT study of the lithium-graphite phase diagram

Graphite is the most widely used and among the most widely-studied anode materials for lithium-ion batteries. With increasing demands on lithium batteries to operate at lower temperatures and higher currents, it is crucial to understand lithium intercalation in graphite due to issues associated with lithium plating. Lithium intercalation into graphite has been extensively studied theoretically using density functional theory (DFT) calculations, complemented by experimental studies through X-ray diffraction, spectroscopy, optical imaging and other techniques. In this work, we present a first principles based model using DFT calculations, employing the BEEF-vdW as the exchange correlation functional, and Ising model to determine the phase transformations and subsequently, the thermodynamic intercalation potential diagram. We explore a configurational phase space of about 1 billion structures by accurately determining the important interactions for the Ising model. The BEEF-vdW exchange correlation functional employed accurately captures a range of interactions including vdW, covalent and ionic interactions. We incorporate phonon contributions at finite temperatures and configurational entropy to get high accuracy in free energy and potentials. We utilize the built-in error estimation capabilities of the BEEF-vdW exchange correlation functional and to develop a methodological framework for determining the uncertainty associated with DFT calculated phase diagrams and intercalation potentials. The framework also determines the confidence of each predicted stable phase. The confidence value of a phase can help us to identify regions of solid solutions and phase transformations accurately.


I. INTRODUCTION
Graphite is the most ubiquitous anode material used in Li-ion batteries owing to its very low potential, high cycle life and low cost. 1,2 Owing to its enormous importance, Li intercalation in graphite has been the subject of numerous experimental and thereotical investigations. [3][4][5][6][7][8][9][10][11] The various phases of Li intercalation into graphite have been characterized using X-ray diffraction (XRD), 12-14 optical methods. 15 However, due to the high degree of disorder, its difficult to precisely pinpoint the structure for a particular phase in the Li-graphite phase diagram, especially the low Li content phases. 3,5 Theoretical studies employing density functional theory calculations have been used to complement the experimental understanding and have painted a clearer picture of Li-graphite phase diagram. 16 However, Persson et al. utilized density functional theory calculations that did not explicitly account for van der Waals (vdW) interactions and used a simple constant vdW correction for the various phases. 16 Graphite comprises of graphene sheets stacked on each other, which are bound by weak van der Waals forces. As Li intercalates in the space between theses sheets, there is an increase in the Li-C covalent interactions along with a decrease in the van der Waals interactions between the graphite sheets. 17,18 Given that the subtle interplay between these two interactions determines the phase diagram, it is crucial to employ an exchange correlation function that is capable of accurately predicting over a wide range of bonding environments. Recently developed exchange correlation functional, BEEF-vdW is designed such that it minimizes prediction error for a range of data sets involving molecular formation energies, molecular reaction energies, molecular reaction barriers, noncovalent interactions, solid state properties, and chemisorption on solid surfaces. 19 In addition, the functional possesses bayesian error estimation, which is designed to reproduce known energetic errors by mapping the uncertainties on the exchange-correlation parameters. This capability allows an error estimation capability and has been employed successfully to understand uncertainties associated with reaction rates in heterogeneous catalysis, 20 activity relationships in electrocatalysis, 21 and mechanical properties of solid electrolytes. 22 We argue that it is crucial to utilize an exchange correlation functional that explicitly accounts for van der Waals interactions and motivate this through a simple example. The adhesion of two graphene sheets is incorrectly predicted to be unstable by 0.13 eV using the PBE exchange correlation functional, while BEEF-vdW correctly predicts that they are stable by 0.07 eV.
In this study, we build on the pioneering work by Persson et al. to provide a refined picture for Li-graphite phase diagram employing density functional theory calculations using BEEF-vdW exchange correlation functional. We apply the cluster expansion formalism to explore the phase space of Li and natural graphite based compounds. The coefficients of the cluster expansion are used to carry out a rigorous search over the enormously large configurational phase space. We explore the large phase space by varying the Li concentration in different unit cells of graphite of multiple sizes. We precisely determine the Li-C, C-C and different Li-Li interactions along with the uncertainties associated with these interaction coefficients. We show that the in-plane Li-Li interaction is dominantly electrostatic in nature. We also incorporate the configurational entropy which plays a significant role in phases with low Li concentration. Based on our analysis, we identify the following stable phases in the Li intercalation diagram viz. Li 0.047 C 6 ,Li 1/6 C 6 ,Li 1/3 C 6 ,Li 2/3 C 6 and Li 3/4 C 6 . Finally, we show for the first time, propogation of uncertainty associated with the density functional theory calculations through the cluster expansion and ultimately to the calculated thermodynamic intercalation diagram. We believe the proposed framework of estimating uncertainty will become widely important in computational investigations for phase diagrams in batteries. We also believe that such an approach utilizing an exchange correlation functional in diverse bonding environments is crucial for understanding anion intercalation in graphite, 23,24 alkali metal (Li, Na, K) intercalation in organic compounds such as melanin, 25 Li intercalation in layered compounds such as MoS 2 , LiV 0.5 Ti 0.5 S 2 , etc. 26,27

II. METHODOLOGY
In this work, we will focus on the electrochemical Li intercalation in natural graphite given by: The phase transformation of Li-graphite compounds is determined by the Gibbs free energy change associated with this process given by: where G LixC6 is the free energy of the Li-graphite phase, G C6 is the free energy of the graphite phase, G Li + is the free energy of the Li ion solvated by the electrolyte and G e − is the free energy of the electron at the potential of the graphite electrode. The sum of the free energy of the Li ions and the electrons can be related to the free energy of bulk Li metal as shown through the reaction: which gives us the relation G Li + + G e − U=0V = G Li (s) . This is termed as the computational lithium electrode and provides a tractable way to determine the sum of the free energies of Li ion and electron for concerted Li ion-electron transfer reactions. 28 Through this relation, the free energy of an electron is now calculated relative to the potential of Li/Li + redox couple, G e − = G e − U=0V − eU Li/Li + . Thus substituting the relation in Eq. 3, we get: Eq. 4 gives the intercalation potential of a particular phase of the Li-graphite phase diagram. To derive the thermodynamic intercalation potential diagram, we need to consider phase transformation from one stable phase to another as Li insets into graphite. The potential for phase transformation from Li x0 C 6 to Li x1 C 6 can thus be expressed as: To calculate the intercalation potentials, we need to calculate the free energies of the stable phases. For a particular phase Li x C 6 , we have large number of possible structures due to a number of possible sites that Li atoms can occupy between the graphene sheets with different free energies. Under a thermodynamic formulation, intercalation will proceed through phases with the minimum free energy. The gibbs free energy comprises of two components namely the enthalpy and the entropy as given by, ∆G = ∆H − T∆S. In the next two sections we will derive the enthalpy and entropy for the possible structures of Li-graphite.

A. Enthalpy Calculation
The enthalpy of formation, ∆H of a system consists of the internal energy change of the system, ∆U and the pressure-volume work, ∆PV. It has been shown that the pressure volume work is negligible compared to the change in internal energy for Li-graphite structures. 29 The internal energy of all Li-graphite structures, Li and graphite is calculated using density functional theory (DFT) calculations. The density functional theory calculations are carried out using GPAW, which is a real space implementation of the projector-augmented wave method. 30,31 The DFT energies are then used to find the enthalpy of formation for each of the structures given by Li intercalation in graphite is an interplay between Li-C interactions, Li-Li repulsion and C-C vdW forces. Thus, it is important to choose an appropriate exchange correlation functional for the DFT calculations. Experiments have shown that Li intercalates in graphite in stages. 3,5 A stage 'n' Li-graphite phase implies that there are n layers of graphene between two adjacent Li layers. In this work, we will explore phases with stages 1, 2, 3 and 4. Phases with higher stages are computationally very expensive to calculate due to the large number of atoms in the unit cell. All DFT calculations were done with BEEF-vdW exchange correlation functional and used a Monkhorst-Pack grid for the Brillouin zone sampling. 32 A convergence of 5 meV/6C was achieved with 10 × 10 × 10 k-point grid for a graphite unit cell (4C atoms). For the Li-graphite phases, the k-points were appropriately scaled down as per the size of the unit cell for each of the phases, so as to maintain the same level of convergence. A grid spacing, h = 0.18Å and a Fermi-Dirac smearing width of 0.05 eV was used.
The graphite AA and AB stacking structures were considered and the AB stacking structure is more stable by 0.09 eV compared to the AA stacking structure, consistent with experiments. 17,18 We also considered the energy for fully intercalated Li-graphite phase (LiC 6 ) for the AA and AB stacking structures and the AA stacking was stable while the AB was unstable as observed by earlier works. 17,18 Thus henceforth, for Li-graphite phases, we assume that the Li intercalated graphene layers will be AA stacked, while the non Li intercalated layers would be AB stacked.
In a recent work, Christensen et al. showed that systematic errors in the formation energy of alkali oxides can be reduced through the use of a reference compound that has the same oxidation state. 33 They demonstrated this through a reduction in the uncertainty as estimated through the BEEF ensemble. In this work, we are interested in investigating Li intercalation into graphite and given the similarity in potentials between graphite and Li metal, we choose Li metal as the reference for all the uncertainty estimates provided in this work. It is worth pointing out that the uncertainty estimates for Ligraphite phases are much lower using Li metal reference than LiCl reference. The uncertainty estimates follow the standard BEEF ensemble estimation procedure and more details regarding this can be found elsewhere. 20,21 As suggested by Christensen et al., we correct for the value of the reference energy of Li through the following equa- LiCl is the experimentally measured formation enthalpy of LiCl.
The insertion of Li into graphite leads to expansion of the lattice. 4 To get accurate lattice constants for each of the Li-graphite structures, we first expand graphite inplane followed by an expansion out of plane and carry out an energy minimization. As mentioned before, we restrict our DFT calculations to calculate the energy of stage 1 to stage 4 structures, since the unit cells become too large for structures beyond stage 4. Further, beyond stage 4 structures have not been clearly observed in experiments. 3 The Li-graphite phases have a large number of structural possibilities due to the large number of sites available for Li in graphite. A systematic way of accounting for this structural disorder is through the method of cluster expansion. Cluster expansion has been previously used to describe the Li-graphite phases and other intercalation compounds like Li x CoO 2 , Li x TiS 2 , etc. 16,34,35 In cluster expansion, the Li sites are considered as a lattice model and an occupation variable is assigned to each Li site. We define the occupation variable, s, to be 1 if the site is occupied by Li and 0 when the site is empty. The formation enthalpy of the system is then derived as a polynomial expansion of these occupation variables as where s i is the occupation variable associated with site 'i', V is the energy associated with Li occupying a site 'i' in the graphite and J i,j are the various interactions between two Li atoms occupying two sites in the lattice model. We assume that three-body and higher order interactions between Li's are negligible. Thus, we limit our cluster expansion to two-body interactions.
The various kinds of two-body interactions J i,j can be associated with the different distances between the Li occupied sites. As these Li-Li interactions are electrostatic in nature, we expect them to decrease with increasing distance. Here, we consider interactions with distance less than 5a where a is the distance between two adjacent carbon atoms in the graphite network. As we show later, interactions beyond this distance are negligible. For structures with stage n ≥ 2 , the distance between occupied sites in different planes is greater than 5a and hence there would be no out-of-plane interactions. However there will be interactions for stage 1 structures. As a result, we have two different cluster expansions, one for stage 1 and another for stage n ≥ 2. One thing of note is that ∆H = 0 for graphite which implies that C = 0 in Eq. 6. Thus the final cluster expansions are given as: In Eq. 7, J ik represent the various Li-Li interactions within a single Li layer (in-plane) as shown in Fig. 1(a). In Eq. 8, V i and V o represent the Li-C interactions for the lower and upper Li layers in the stage 1 unit cells respectively, while the J ik and J ok represent the in-plane and out of plane interactions as shown in Fig. 1(b). We count all interactions in the unit cell and with neighboring unit cells with appropriate weights. For the stage 3 and stage 4 structures, out-of-plane interactions are negligible and hence neglected, which means that their energy can be calculated from the stage 2 structures as follows: (10) The last two terms in Eq. 9 and Eq. 10 account for the energy required to reorient the grahene sheets to a different stacking. To determine the coefficients of the cluster expansion, we calculated 32 structures using DFT such that atleast one structure had each of the interactions and all structures were unique. A regression analysis on these 32 DFT calculated structures gives all the interactions. Using the values for cluster expansion coefficients, the formation enthalpy for the entire phase space can be determined. We explore the phase space through various unit cell sizes of the graphite where we intercalate Li atoms and calculate energies through cluster expansion. We could explore unit cells with a maximum of 24 Li sites for stage 1 and stage 2, beyond which the combinations of filling Li atoms become very large (∼ 2 24 ). Thus we can only explore upto the phase Li 0.047 C 6 . Ordered Li x C 6 phases with x < 0.047 have not been observed experimentally due to the formation of SEI during the initial electochemical Li intercalation. 5   FIG. 2. The increase in the C-C interlayer distance among adjacent graphene sheets as a functions of the various stages of Li intercalation. Note that the interlayer distance increases non linearly proving that the vdW interactions get screened non uniformly upon intercalation.

B. Entropy Calculation
The entropy change associated with the formation of the Li-graphite phase can be written as ∆S = S LixC6 − S C6 − xS Li . Entropy of a material consists of the configurational entropy and the entropy of the vibrational modes of the structure. Here we assume that the entropy of vibrational modes is approximately the same as the entropy of vibrational modes of graphite. Similar conclusions have been reached for other vdW bonded structures upon intercalation. 25 Thus the entropy change is now given by: The configurational entropy of Li x C 6 will be derived later in the results section. The standard state entropy of Li is 29.09 J/mol K. 36 The formation enthalpy and entropy are then used to find the Gibbs free energy for all phases. The free energy diagram is used to identify the most stable phases and the phase transition regions which in turn give the intercalation potential diagram for Li intercalation in graphite.

III. RESULTS
The formation enthalpies for 32 structures have been calculated using DFT as per the procedure described in the Methods section. Table II in the Appendix provides the different kind of interactions, the exact Li x C 6 phase, the lattice constants for each for each of the structures, the stage number and the formation enthalpy for the 32 structures considered.
The in-plane lattice constant a increases slightly upon Li intercalation. However, out-of-plane lattice constant, c, which is the distance between adjacent graphene sheets  significantly increases upon Li intercalation to a maximum increase of 11.16% upon full intercalation. It is worth pointing out that c is more strongly dependent on the stage number of Li x C 6 than the filling fraction, x. A plot of c as a function of the stage is shown in Fig. 2. This implies that the upon intercalation, the Li is weakening the vdW interactions between the graphene planes. It is also worth highlighting that this plot shows that the weakening effect is non-linear and hence a simple constant vdW correction as assumed by Persson et al. 16 is prone to errors. Using the kind of interactions involved in each structure along with its formation energy, we determined the interaction coefficients for the cluster expansions for stage 1 and stage n ≥ 2. We use the least squares optimization method to evaluate the interaction coefficients. The maximum error for the fit for stage 1 cluster expansion is 0.047 eV and for stage n ≥ 2 cluster expansion is 0.05 eV. These error are maximum deviations of the cluster expansion model from the DFT energies after performing the least squares method analysis to determine cluster expansion coefficients using MATLAB software. The evaluated interaction coefficients are shown in Table I.
As we can observe from the table, the occupation terms, V, which correspond to Li-C interactions are negative implying that these interactions are attractive while the J ii s and the J oi s are positive implying that the Li-Li interactions are repulsive in nature. The in-plane Li-Li interactions (J ii ) are plotted against the inverse of the corresponding distances of the interactions and shown in Fig. 3. The in-plane Li-Li interactions decrease with increasing distance and are proportional to 1/r showing that the interactions are dominantly electrostatic in nature. However, as can be seen from Table, the out-ofplane Li-Li interactions (J oi ) are small and do not follow a similar trend with distance associated with interactions. This is due to varying degrees of charge screening by the graphene sheets for these interactions.
We use the procedure described in the methods section to calculate the formation enthaply for different unit cell sizes upto 24 Li sites for different phases of Li x C 6 and different stages. We calculated formation enthalpies for about 26000 different structures. To determine the free energy of structures we also need to calculate the entropy of these structures. As shown in the methods section, the entropy depends on the phase Li x C 6 and the configurational entropy associated with that particular phase.
To derive the configurational entropy, we can consider a graphite block with N atoms. A graphite block of N atoms due to its hexagonal symmetry will have N/2 sites for Li to occupy. From the interaction coefficients, we can see that the first in-plane interaction coefficient is significantly larger than the others, which implies that it is very unfavorable for two Li atoms to occupy adjacent sites. Hence due to the hexagonal symmetry of graphite, we can only fill a Li among 3 adjacent sites. This implies that for a particular phase x of Li x C 6 , the total configurations Ω of possible Li ordering is the number of ways of arranging Nx 6 Li atoms in N 6 sites. Thus Ω = N 6 Nx 6 . The configurational entropy for Li x C 6 for N carbon atoms can now be written as follows: For our model, we need to normalize the configurational entropy to 6 carbon atoms. i. e. s conf = 6 N S conf . Now we can derive the expression for configurational entropy of phase x as given by: After applying the Stirling's approximation for large N, the configurational entropy can be expressed as: The formation entropy then can be evaluated for each of the structures as From Eq. 15, we evaluate the entropic contribution to free energy formation of Li x C 6 at T = 298 K comprising of the configurational and the entropy change of the Li atoms in graphite relative to lithium metal. This entropic contribution for various phases x of Li x C 6 is shown in Fig. 4 and it agrees well with experimental measurements done by Reynier et al. 5 In the work by Reynier et al., the entropy change was calculated by measuring the change in open circuit voltages in half cell cells with graphite electrodes for a small change in temperature.
The enthalpy and entropy at standard temperature of 298 K were combined to determine the Gibbs free energy for 26000 different structures as shown in Fig. 5(a). To construct the phase diagram of Li-graphite, we need to choose the structures which lie on the convex hull of free energy vs. phase x diagram. The phase diagram being the convex hull is a result of the condition that a system in thermodynamic equilibrium minimizes its free energy. The convex hull is derived from the 26000 points which gives the phase diagram as shown in Fig. 5(b). At all the stable phases, the structures within an energy of k B T were considered equivalent. From the convex hull, we evaluate the intercalation potentials, using the free energy of the most stable phases. The intercalation potential diagram is shown in Fig. 6. The calculated intercalation potentials are within 0.2 V of the experimentally measured potentials also shown in Fig. 6. It is worth pointing out that the comparison between experiment and theory is within the uncertainty associated with the DFT-calculated intercalation potentials. The uncertainties for the free energy of formation were evaluated by calculating the free energies for an ensemble of 2000 exchange correlational functionals and then taking the standard deviation. The uncertainty for the intercalation potentials was evaluated by calculating the free energies of the stable phases and in turn the intercalation potentials with 2000 exchange correlational functionals and then determining the standard deviation.
From our phase diagram, we find that at low Li concentration phases, Li x C 6 is a mixture of all 4 stages with stage 4 being most stable for x ≤ 1/6, then it phase transitions to a mixture of stage 3, 2 and 1 for 1/6 < x ≤ 1/3, followed by a mixture of stage 2 and 1 for 1/3 < x ≤ 1/2 and then finally to a fully stage 1 compound for x > 1/2. All phase transitions and the interacalation potentials associated with these phase transitions agree very well with those observed in experiments by X-ray diffraction, 3 spectroscopic 14 and optical measurements 15 especially for x ≤ 0.5. The vdW forces and the Li-Li repulsive interactions and their non linear variation with the x (of Li x C 6 ) play a key role in determining the phase transitions. We also report the uncertainties associated with the free energies for each of the stable phases as calculated from DFT. Theses uncertainties are generated through calculations of the energy through the generation of an ensemble of functionals as implemented in BEEF-vdW. 19 The uncertainties have been used to determine the uncertainty associated with the DFT-calculated thermdynamic intercalation diagram. Finally, we would like to re-emphasize the importance of vdW interactions and the need to appropriately choose the exchange correlational functional to determine accurate energetics.
The theoretical results from our approach agree very well with experiments in terms of the structures and energies of the stable phases as well as the intercalation potential diagram. We have clearly demonstrated that the subtle vdW interactions need to be captured accurately by choosing an appropriate exchange correlational functional in DFT as they determine the staging and the free energies of the phases, which then determine the intercalation potential diagram.

IV. SUMMARY AND CONCLUSION
The presented framework comprising of DFT and statistical thermodynamics accurately predicts the phase diagram and intercalation potentials for electrochemical Li intercalation in natural graphite. We highlight the importance of accurately capturing vdW interactions in the exchange correlational functional employed in DFT and we show that this is crucial to accurately determine the stable phases during intercalation. We observe stable phases at x = 0.047, 0.167, 0.333, 0.5, 0.667, 0.75, 1. We also observe the change in staging from stage 4 to stage 1 for these phases with increasing Li inter-calation, in good agreement with experiments. Cluster expansion technique and DFT calculations with the correct exchange correlational functional play a key role in precisely determining the C-C, C-Li and Li-Li interactions. We show that the in-plane Li-Li interactions are predominantly electrostatic in nature. Finally, we introduce a method to evaluate uncertainty of free energies of various phases and the associated intercalation potentials through the built-in error estimation capability of the BEEF-vdW exchange correlation functional. Lastly, we would like to point out that vdW and entropic contributions, even though small, play a key role for determining intercalation processes in various anode materials where the intercalation potentials are close to that of Li/Li + .