On-the-fly assessment of diffusion barriers of disordered transition metal oxyfluorides using local descriptors

Disorder plays an increasingly important role in the design and development of high-performance battery materials and other clean energy materials like thermoelectrics and catalysts. However, conventional computational design approaches based on the thermodynamic properties of statistically averaged structures are unable to predict the accessible energy and power densities of such materials. Kinetic properties like ionic diffusion within locally resolved atomic structures is needed to perform longer time and length scale simulations like kinetic Monte Carlo in order to accurately estimate kinetic properties like power densities in battery electrodes. Here, we present and demonstrate a fast, on-the-fly, approach to calculate local diffusion barrier as a function of only the local atomic structure using machine learning and cluster expansion, particularly for Li-ions in lithium-rich transition metal oxyfluorides and the disordered rock salt (DRS) Li2–xVO2F electrodes.

• Lithium diffusion in disordered transition metal oxyfluoride battery electrodes • Fast on-the-fly prediction of diffusion barriers with machine learning • Simple local structural descriptors for kinetically resolved activation barriers • Discovering kinetic barrier descriptors with LASSO

Introduction
Lithium-rich transition metal oxyfluorides have garnered significant interest in recent years for their application as high energy density cathode materials for lithium-ion batteries. These materials are typically synthesized using a high-energy ball-milling approach, which introduces a high degree of site disorder [1,2]. Lithium-rich transition metal oxyfluorides have a disordered rock salt (DRS) structure where lithium and transition metal occupy cationic sublattice, while oxygen and fluorine occupy the anionic sublattice. Many studies of DRS oxyfluorides have followed the original experimental report on the promising performance of Li 2 VO 2 F in 2015 [3]. The primary research focus has been on improving its cyclability while retaining its high energy density. The impact of substituting vanadium with different transition metals (e.g., Co, Fe, Mn, Mo, Nb, Ti, W), the extent of over-lithiation and oxygen and fluorine ratio have been investigated [2,4,5,6,7,8]. Despite the continued effort, the inner workings of DRS oxyfluoride remain largely elusive, necessitating further research to gain a sufficient understanding of their thermodynamic and kinetic properties. The insights gained on these disordered materials will have impact that reach far beyond battery research, as disorder is present in materials for other clean energy applications such as thermoelectrics and catalysts.
The thermodynamic properties of the DRS transition metal oxyfluorides can be probed computationally. A commonly used approach is to use Cluster Expansion (CE) in conjunction with first-principles simulations such as density functional theory (DFT) calculations [9]. The computational investigations based on CE were successful in predicting crucial aspects of DRS materials such as order-disorder phase transition and open-circuit voltage profile [10,11]. While providing crucial insights, these models largely overlook the kinetic contributions (i.e., transport of lithium ion upon charge and discharge). Consequently, the model predictions are valid at the thermodynamic limit, where the battery is charged/discharged at a very slow rate. This limitation poses a challenge in evaluating material's accessible energy and power density as they are largely governed by the kinetic properties of the material. Furthermore, the kinetic properties are expected to play a pivotal role in determining the cyclability of the material.
The transport kinetics of cathode materials are often assessed using the energy barrier of the diffusing lithium ion when moving to an adjacent site. Lithium ions in DRS transition metal oxyfluorides predominantly move from an octahedrally coordinated site to another as shown in Fig. 1. Two adjacent octahedral sites are joined by two tetrahedrons. Lithium diffusion in DRS occurs either in a straight line between two anions or through one of the two tetrahedrons. The transport in a straight line is referred to as an oxygen dumbbell hop (ODH) as the term was originally coined for oxides, and the transport through the tetrahedron is referred to as tetrahedron site hop (TSH) [12,13]. The transport mechanism is largely governed by the local atomic surroundings of the diffusing lithium ion. One of the key characteristics that favor macroscopic lithium diffusion is the existence of a percolation path formed by consecutive TSH with no transition metals in its surrounding octahedra (so-called 0-TM path). This condition can be achieved by having a high degree of fluorination (above 15 %) [7], as the amount of fluorination should be matched by the over-lithiation to keep the charge neutrality of the system. By contrast, a low degree of fluorination does not allow the formation of the 0-TM path, and its impact is largely dominated by its local chemical compositions of TSH. Indeed, it has been found that the low degree of fluorination (below 15 %) can have a negative impact on the lithium diffusion [7]. Figure 1: Diffusion mechanism of octahedrally coordinated lithium ion. Yellow sphere represents diffusing lithium ion while red and blue spheres represent anions (oxygen or fluorine) and gate cations (lithium, transition metal or vacancy). The vacant lithium-ion site is shown as a circle with dashed lines.

TSH
Euchner et al. have recently investigated the impact of the local surroundings of the diffusing lithium ion in both DRS oxides and oxyfluorides [14]. The diffusion mechanism description is simple for layered oxides; i.e., the ODH mechanism prevails for mono-vacancy diffusion (none of the gate-cation sites are vacant), while the HSP mechanism prevails for di-vacancy diffusion (one of the gate-cation sites is vacant). In addition, the di-vacancy diffusion has significantly lower diffusion barrier as it allows for more open space and reduces electrostatic repulsion. However, the simple diffusion mechanism description is not applicable for DRS oxides and oxyfluorides, and the correlation between the diffusion barrier and the local atomic arrangement is not straightforward. Nevertheless, the thorough study of Euchner et al. [14] revealed that (1) substituting gate-cation sites occupied by a transition metal with a lithium atom lowers the diffusion barrier, (2) substituting oxygen atoms in dumbbell sites (two anions in the middle of the straight diffusion path) with fluorine increases the diffusion barrier, (3) substituting other ligand oxygen atoms could increase or decrease the barrier depending on their sites. It was also found that the disorder shifts the TSH mechanism to become energetically favorable compared to the ODH mechanism even for mono-vacancy diffusion.
There are still several open questions that remain to be answered: How do we determine if the diffusing lithium atom will pass through upper or lower tetrahedron? How do non-local quantities such as the lithiation level influence the diffusion barrier? Can we develop a simple and fast computational framework to estimate diffusion barriers using descriptors that are easy to quantify? The main objective of this work is provide answers to the question above, thereby allowing one to incorporate on-the-fly estimation of lithiumion diffusion barrier in DRS transition metal oxyfluorides, e.g., to optimize and extend regions of high power operation. Such quick barrier estimation will pave the way for developing an efficient kinetic Monte Carlo simulation scheme that accounts for kinetic contributions when simulating the structure evolution. In this work, we report our work on estimating the diffusion path (up, down or straight) and the diffusion barrier using regression model fitted to the results obtained using the nudged elastic band (NEB) method based on DFT calculations. We examine Li 2 -x VO 2 F as an exemplary case, and our results agree well with previously reported observations and show which descriptors play a critical role in determining the diffusion barrier.

Theory
All of the diffusion barriers used for fitting our models are calculated using DFT performed with the Vienna Ab initio Simulation Package (VASP) [15,16,17,18]. The projector augmented-wave (PAW) method was used to describe electron-ion interaction [19]. All of the DRS Li 2 -x VO 2 F were created using 3 × 4 × 4 supercell containing 96 atoms in case of no vacancies. The generalized gradient approximation as parametrized by Perdew, Burke and Ernzerhof [20] was used as the exchange-correlation functional, while a rotationally invariant Hubbard U correction [21,22] was applied to the d orbital of V with the U value of 3.25 eV. The plane-wave cutoff of 500 eV was used. Both cell and atomic positions were fully relaxed such that all the forces are smaller than 0.02 eV/ A for the initial image of the NEB, while the same settings were used for the other images except the cell shape was fixed to that of the initial image. Integration over the Brillouin zone was carried out using the Monkhorst-Pack scheme [23] with a 3 × 3 × 3 k-point grid. A total of 48 climbing image NEB simulations were performed, while each simulation had 5 intermediate images between the initial and final images. Access to the dataset is provided through the Materials Cloud Archive (insert link prior to publication).

Features for the lithium diffusion model
A total of 11 features are devised and investigated to fit the lithium diffusion characteristics calculated using the NEB method. The geometric features are shown in Fig. 2. Consequently, we employ a total of 9 local, geometry-based features (5 for anions and 4 for cations). Diffusion characteristics are also affected by some non-local contributions [14], and two non-local features are also introduced. The non-local descriptors are: the energy difference between the initial and final states of diffusion (∆E = E final − E initial ) and the lithiation level of the compound. It is noted that one can extract the values of all the features quickly except for ∆E, which is typically determined using time-consuming DFT calculations. However, the ∆E value can also be estimated quickly with a good accuracy using CE, which will most likely be used in conjunction with the diffusion barrier estimation scheme presented in this work. The description of all features employed in this work are summarized in Table 1.

Estimation of diffusion path
The first crucial step for estimating the diffusion characteristics is to determine which diffusion path it takes: TSH through the top tetrahedron, TSH through the bottom tetrahedron or ODH. The distinction of the path The linear regression coefficients of the 3 features are shown in Fig. 3a. It can be seen from the coefficient values that the vacancy and F atoms "attract" the lithium atom in the transition path, while the V atoms "repel" the lithium atom. The vacancy has the greatest impact by far, which is in good agreement with the well-known fact that the vacancy is the main cause of the TSH mechanism and significantly lowers the energy barrier [12,13]. The V substitution of Li in the gate cation sites is also reported to increase the energy barrier due to its increased electrostatic repulsion [14]. Fluorination Lithium concentration, as ratio between number of Li and number of F atoms the whole system. is reported to have a small negative impact on lithium diffusion [6], which also agrees with the coefficient value. The predicted diffusion heights are plotted against the heights from the NEB calculations as shown in Fig. 3b. As expected, the predictions of the exact height is not very accurate with the cross-validation RMSE of 0.64 A. However, the objective of fitting the height of the path is simply to determine which tetrahedron it passes through. In that regard, the model predicts the correct path with 87.5 % accuracy, while the model predicts the path to go through a wrong tetrahedron or through the dumbbell anions in 1/8 of the structures. The 3 predictions that are exactly 0 are caused by the nonuniqueness of the features, i.e., in these cases the top and bottom paths are equivalent with respect to the features. The overall success rate of 87.5 % is satisfactory for such a simple linear model, and it can be used to classify   which cations and anions play a primary role in determining the energy barrier.

Estimation of diffusion energy barrier
Instead of targeting the forwards (from the initial configuration to final configuration) and backwards (from the final configuration to the initial configuration) energy barriers, a common practice is to aim for the so-called kinetically resolved activation (KRA [13]) barrier, which corresponds to the energy difference between the transition state energy and the average energy of the initial and final configurations (see Fig. 4). In this work, we follow both approaches, building first a model that aims to calculate the forward barrier, and later on a second model that targets KRA.
We fit linear regression models to the forward energy barrier using the 11 features in Table 1 and the whole dataset of 48 barriers. We note that the diffusion paths are oriented such that the diffusing lithium passes through the top tetrahedron. The coefficients of the linear fit are shown in Fig. 5a. To assess the prediction accuracy of the models, we perform leave-one-out cross-validation, i.e., for every example in the dataset, we fit the 11 features to the 47 other data examples and predict the value of the single example. The predicted barriers are plotted against the NEB-calculated barriers in Fig. 5b. The first model aiming for the forward barrier ( Fig. 5a and Fig. 5b) predicts a 0.44 coefficient for the ∆E, which is consistent with the KRA approach where the weight of ∆E is fixed to be 1/2. The close agreement

Energy
Reaction coordinate initial final TS of the weight of ∆E using our approach with that of KRA demonstrates its robustness, despite its simplicity.
We now move on to the fitting results using the KRA model and determine which features are most relevant for describing the energy barriers. We fit another set of linear models with a 1 penalty on the weights to determine the most relevant features to include. These models are known as least absolute shrinkage and selection operator (or LASSO) [24]. The LASSO cost function is: where y is a R N vector with target values (energy barriers in our case), X ∈ R N ×F is a matrix where each row corresponds to one training example and each column corresponds to the features in Table 1. The objective is minimized with respect to the weight vector w ∈ R F and α is a hyperparameter, which controls the trade-off between optimizing the model fit and minimizing magnitude of the weights. A high α value will force many of the w entries to 0 and thus induce sparsity. The features are standardized before fitting the KRA model, i.e., rescaled and shifted to zero mean and unit variance. This allows for a direct comparison between the corresponding linear regression coefficients. We minimize the LASSO cost function (1) on the whole dataset (N = 48) and use 20-fold cross-validation to determine the optimal value of α. The KRA model is implemented using the LassoCV class in scikit-learn. The coefficients of the KRA model is shown in Fig. 5c. The coefficients of ∆E, the gate anion and initial and final ligands are zero, indicating that they have zero contribution To simplify the model we take the relevant features determined by the LASSO method and refit a linear model using only these five features. The result is shown in Fig. 5e and Fig. 5f. In this case we do not standardize the features, and therefore the model coefficients are directly interpretable, e.g., the coefficient corresponding to GCT X shows that by replacing one of the top gate cations with a vacancy, decreases the barrier with approximately 0.3 eV and replacing O with F in top of the tetrahedron increases the barrier by approximately 0.06 eV. The final CV RMSE is 0.12 eV. Despite its simplicity, the average error of our final model, 0.12 eV, is in a similar range to the energy prediction error one would expect from a CE model for the initial and final configurations. Therefore, the model can be used in conjunction with the CE approach to develop an on-the-fly kinetic Monte Carlo model to assess the formation of percolating network and to evaluating the accessible energy and power density of the DRS transition metal oxyfluorides.

Conclusion
This work provides a pathway to fast, on-the-fly, estimation of the power density of disordered battery electrode materials, where local diffusion kinetics control the bulk properties instead of the average thermodynamics of different structural motifs. We show, that a combination of, e.g., cluster expansion (CE) for energy evaluation and our ML-based kinetic barrier prediction model would allow access to long time/length scale simulations towards quickly assessing the power density of the materials taken disorder in to account.
We have defined a number of local features, from which the barriers are estimated using linear models. The most relevant features are identified using LASSO regression and used in a simplified model. Our on-the-fly model has an RMSE of 0.12 eV for kinetically averaged barriers. The uncertainty in the full ionic diffusion kinetics is also contributed from the uncertainties in energies predicted from CE (or DFT for that matter), which are typically on the order of 0.1 eV [11].
The experimentally/practically observed lithium diffusion properties at longer length scales originate from percolation of lithium diffusion paths that are just large enough (i.e., the percolation threshold) to create a diffusive pathway throughout the macroscopic system [25]. We found that the presence of vanadium as a gate cation has a negative impact on lithium diffusivity, which aligns well with the 0-TM percolation network model [7]. We also found that the presence of fluorine ions in the tetrahedron (TAT F) has a slightly negative impact on lithium diffusivity, which agrees with the experimental observations of the hindrance of lithium diffusion at low fluorination levels, where the 0-TM percolation network is not yet formed [7]. Overall, fluorination can facilitate the lithium diffusion on a macroscopic scale by forming the 0-TM percolation network, while still having a slightly negative impact on a local, site-to-site hopping level.
The nature of uncertainties from both the ML models and the CE method are such that those are larger when the barriers themselves are larger. This further ramifies the impact of our approach as a valuable design tool, as the highest performance materials are expected to have small percolation threshold for the diffusion barrier and we need to be accurate on only the barriers smaller than that to have a good description of the system. In short, our approach is reliable where it matters the most, such that we can screen systems quickly, reliably and definitively for power capabilities.

Acknowledgments
Authors acknowledge support from the European Union's Horizon 2020 research and innovation program FET-OPEN project LiRichFCC under Grant Agreement #711792 and the BIG-MAP project under Grant Agreement #957189, as well as the Villum Foundation's Young Investigator Programme (4th round, project: In silico design of efficient materials for next generation batteries. Grant number: 10096). We hereby submit the revised version of our manuscript entitled "On-the-fly assessment of diffusion barriers of disordered transition metal oxyfluorides using local descriptors" for your consideration for publication in the ISE-2020 special issue of Electrochima Acta.