Bayesian-Driven First-Principles Calculations for Accelerating Exploration of Fast Ion Conductors for Rechargeable Battery Application

Safe and robust batteries are urgently requested today for power sources of electric vehicles. Thus, a growing interest has been noted for fabricating those with solid electrolytes. Materials search by density functional theory (DFT) methods offers great promise for finding new solid electrolytes but the evaluation is known to be computationally expensive, particularly on ion migration property. In this work, we proposed a Bayesian-optimization-driven DFT-based approach to efficiently screen for compounds with low ion migration energies (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\boldsymbol{E}}}_{{\boldsymbol{b}}}{\boldsymbol{)}}$$\end{document}Eb). We demonstrated this on 318 tavorite-type Li- and Na-containing compounds. We found that the scheme only requires ~30% of the total DFT-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\boldsymbol{E}}}_{{\boldsymbol{b}}}$$\end{document}Eb evaluations on the average to recover the optimal compound ~90% of the time. Its recovery performance for desired compounds in the tavorite search space is ~2× more than random search (i.e., for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\boldsymbol{E}}}_{{\boldsymbol{b}}}$$\end{document}Eb < 0.3 eV). Our approach offers a promising way for addressing computational bottlenecks in large-scale material screening for fast ionic conductors.


Results
Chemical search space and crystal structure description. Tavorite-type compounds with a general formula AMXO 4 Z (A: Li, Na; M: group 2, 3, 4, 13 elements; X: group 14, 15, 16 elements; Z: F, Cl, Br, I) were targeted for the E b -based solid electrolyte screening. The model crystal structure (P 1) is shown in Fig. 1a with the host framework comprising with MO 4 Z 2 octahedra (M 1a, 1h; O 2i) and XO 4 tetrahedra (X 2i; O 2i). The MO 4 Z 2 octahedra are corner-linked together at their trans-Z atoms to form chains along [111]. Each oxygen atoms from these chains are in turn shared with X atoms which then assumes a tetrahedral environment. Site splitting occurs for the A atoms (2i). Overall, the search space includes LiMXO 4 F dataset taken from our previous work 13 and newly calculated datasets for LiMXO 4 (Cl/Br/I) and NaMXO 4 (F/Cl/Br/I). Although there are different local A cation pathways in the tavorite structure, our previous calculations determined that its ionic conduction is anisotropic, with the dominant transport pathway being facile in one major cell direction 13 . This conduction pass is defined by a series of local site-to-site jump environments, each sandwiched between two MO 4 Z 2 octahedra. Hence, E b sampling by NEB method was carried out only at the characteristic local pathway bottleneck, as shown in Fig. 1b in asterisk. Figure 2 shows the distribution of Li and Na DFT-E b datasets (a total of 318 samples) that were prepared in advance for the BO-driven search. We note that although the dataset size may not be large enough for practical material discovery, it should be sufficient enough (considering the heavy computation cost of DFT approach for kinetics-related properties) for evaluating how fast BO-driven search can find the best or nearly-best one in a quinary system. Differences in the sample distribution of the two datasets are revealed by estimating their sample statistics such as maximum

DFT-E b dataset.
, and skewness (α 3 ) and kurtosis (α 4 ): quantities are 594} for Li and Na, respectively. This comparison clearly shows that the Na case has a broader range and more samples with large E b which could mainly stem from the larger atomic mass and ionic radius of Na ( + r Na = 1.02 Å vs. + r Li = 0.76 Å for an octahedral environment) 44 . On another note, both distributions are positively skewed (α 3 > 0) but with the Na case having a heavier tail towards large E b values (α Na 4, > α Li 4, ). These datasets should provide a more stringent performance check for the BO-driven search since random search can favorably sample in the low-E b density region. The optimal compounds ( * x ) in the Li and Na datasets are LiScSbO 4 I (E b = 0.104 eV) and NaErAsO 4 Cl (E b = 0.116 eV), respectively; both are still unreported compounds. Figure 3 shows the schematic workflow for the E b -based BO-driven search within the AMXO 4 Z tavorite search space. At first, the search space of compounds is populated by various ionic substitutions at the A, M, X, and Z sites. Next, t = 5 initial randomly picked compounds ( The predicted favorable conduction pass for A cations within the tavorite framework (A atoms removed) as shown in its 1 × 2 × 2 supercell (in black, along c-direction) (13). The local barrier height E b marked by asterisks are equivalent characteristic path bottlenecks. The VESTA software was uses for structure visualization 43 . which then defines the acquisition function x a( ). Maximizing x a( ) then enables for deciding the next query compound x t to be evaluated for DFT-E b . The sequence is continued until a user-defined number of evaluations or stopping criterion is achieved. In this work, the number of function evaluations was set equal to the number of test data samples.

BO-driven DFT-E b search workflow.
Performance evaluation of BO approach. In this paper, the additive BO model is labeled as aBO while the ordinary BO model is labeled as oBO. Figure 4a shows the efficiency of the three search methods for minimizing the residual gap at each evaluation step t between optimal * x f ( ) and the current best solution ( ) for the Li test data. In the high uncertainty regime (t < 20) of the simulated search (i.e., high σ 2 since majority of test data compounds are still unobserved), aBO shows the best performance. Meanwhile, oBO performs slightly poorer than random search but when t > 20, it starts to outperform. This behavior for oBO especially at the early stage of the search can be explained by its kernel complexity and the skewed E b distribution (see Fig. 3). It should be emphasized though that the nature of the distribution for x f ( ) is usually not known in advance but incidentally  ). The vertical axis in (b) shows the average probability ratio of discovering the optimal Litavorite compound * x . Note that both additive BO (aBO) and ordinary BO (oBO) search methods here used the tuned hyperparameters from half of the Li dataset excluded for search performance comparison (gray area).
Scientific REPoRTS | (2018) 8:5845 | DOI:10.1038/s41598-018-23852-y even with a non-normal distribution the BO approach is still generally efficient in querying for low-E b compounds. Figure 4b shows an alternative performance comparison analysis which is based on the probability ratio of discovering the optimal Li-tavorite compound * x as the number of observations t increases (again, as averaged over 1000 trials). The plot indicates that at t = 35, there is already ~80% probability of discovering * x for both BO-driven searches. Figure 5a shows the search performance of the 3 methods, for the Na dataset: random search, oBO search, aBO search. Two transfer settings for BO were used: transferred hyperparameters only (Li-hp) and with both transferred hyperparameters and posterior GP (Li-GP) from Li dataset for BO. The plots are averaged over 50 trials for a search space of 154-compounds test data. For t < 40, oBO is overall performing poorer than random search and aBO regardless of the inherited model settings. This can be primarily explained in a similar fashion as with the Li case, that is, from the viewpoint of kernel complexity and high estimation error when the number of unobserved compounds is still high. Out of the 5 tested models, Li-GP aBO gains a clear advantage over random search for t > 20. These results validate the use of model transfer and demonstrate the predictive power of the trained GP model from the Li dataset for the Na dataset. Additionally shown in Fig. 5b is the probability ratio of finding the optimal Na-tavorite compound * x . At t = 50 evaluation steps (i.e., ~30% of the search space observed), Li-GP aBO and Li-GP oBO can find * x ~90% and ~80% of the time, respectively. Na and Li ionic conductivity (or Li and Na ion migration energy) are inherently different properties and normally cannot be optimized simultaneously. However, with the systematic approach of knowledge transfer such as in the problem setting (i.e., from Li to Na system), we demonstrated that indeed we can efficiently optimize and find the optimal compound(s) better than random method.
The goal of the BO-driven search can be modified so that compounds that satisfy a cutoff value are explicitly searched, in contrast with just finding the single most optimal compound * x . To demonstrate this, we used the Li-GP transfer model setting and set a criterion of E b < 0.3 eV, referred from perovskite Li 0.34 La 0.51 TiO 2.94 solid electrolyte 45 . The Na dataset was used for performance check and results are displayed in Fig. 6. The vertical axis represents the average number of desired compounds found, which for the Na dataset, would be 17 total compounds meeting the cutoff. For t < 50 (~30% search space coverage), aBO found twice the number of desired compounds than oBO and Random search, discovering ~73% (12.40 compounds) as compared to ~37% (6.26 compounds) and ~33% (5.54 compounds), respectively. However, we note here that aBO failed to find the remaining compounds with E b < 0.3 eV even up to ~80% search space coverage (t = 130). This issue is due to the method trading off some of its predictive accuracy for model flexibility. Still, aBO demonstrates its remarkable performance and suitability for large-scale material screening tasks, given that the search is prioritized on maximizing search hits for desired compounds with as few number of DFT-E b calls as possible.
Descriptor contribution towards E b prediction. Another advantage of additive Bayesian optimization is that the importance of each group of descriptors can be easily interpreted. Figure 7 shows the contributions of descriptor groupings for Li-GP aBO towards E b prediction. The degree of contribution was calculated by taking the normalized ratio of the covariance amplitude σ f 2 for each groups. The two main contributions came from descriptor groups related to the RDF features (g5) and lattice cell features (g1). Meanwhile, inter-polyhedron features (g4) does not contribute and thus could be removed, reducing model complexity from M = 5 down to M = 4 terms. This non-contribution of inter-polyhedron features may be explained by their redundancy since the interatomic-based information contained in them could have been well-expressed already or have been better expressed by RDF features (g5). RDF features, on the other hand, are determined here as effective descriptors for the prediction of E b with an inherently structure-independent nature, making it directly applicable for material search/screening tasks with multiple structure types. ). The vertical axis in (b) denotes the percentage ratio of discovering the optimal Na-tavorite compound * x .
To investigate the characteristics of the descriptor values among different compounds, we analyzed the data distribution of some descriptors. We used g1 descriptors (lattice cell features) which were determined automatically by the present BO method as significantly contributing towards E b prediction (see Supplementary Fig. 4). We observed that for g1 descriptors, there are different distribution shapes, modalities (unimodal, multi-modal) and degree of skewness for the distribution of values which are indicative of variability and variety in the captured information. In addition, the ranges of each descriptor distributions (see Supplementary Table 5) indicate a varying degree of closeness of values among compounds. Nevertheless, g1 descriptors may be generally regarded as sufficiently differentiating for tavorite compounds. As an example, we examined descriptor df which represents the path bottleneck size for the migrating ion. A value of 0.707 Å (minimum among compounds) would make it geometrically unfavorable for Li ion to pass through (Li ionic radius is 0.76 Å in octahedral coordination, as in the tavorite structure), and much more unfavorable for Na ion (1.02 Å) 44 . Meanwhile, a value of 2.240 Å (maximum among compounds) means both Li and Na ion can pass through geometrically.
Based from above importance analysis on descriptor group contributions, we have shown that our chosen set of descriptors and the strategy of grouping them in their natural groups to define sub-kernel spaces for the BO method is indeed an effective approach for navigating the ion migration energy landscape of the tavorite AMXO 4 Z search space.

Post-processing of E b -screened tavorite compounds.
In an actual material screening task, compounds of interest are usually not only evaluated against a single property but also against other metrics. For example, screened compounds after simulated BO can be further narrowed down by thermodynamic stability criterion to assess whether they can be synthesized by experiment or not. For this purpose, we carried out DFT   46,47 . Briefly, the thermodynamic stability energy (E d ) of a given compound was checked against all possible linear combinations of competing phases found in the Materials Project (MP) database 9 . A compound phase may then fall under three cases: (i) E d = 0, the compound is predicted to be at the thermodynamic ground state, (ii) E d > 0, there is a driving force for decomposition, and (iii) for E d ≈ 0, a compound is metastable and may be stabilized by appropriate synthesis condition or high kinetic barriers 5 . Based from this classification and from previous empirical results for DFT formation energies, a value of 0.1 eV/atom was chosen as a reasonable upper limit for stability and metastability 5,8 .  Table S4. Figure 8b,c show the total density of states of LiZrGeO 4 F and NaHfSiO 4 F, with DFT-PBE electronic band gap energies determined to be 4.177 and 4.876 eV, respectively. These values are comparable with other known candidate solid electrolytes such as garnet Li 7 La 3 Zr 2 O 12 (5.79 eV by HSE06) and Li 10 GeP 2 S 12 (3.6 eV by PBE) which have wide band gap, indicative of being able to meet the requirement for very low electronic conductivity 48,49 . Additional data are provided in Table S3 for DFT-optimized structural information. The ionothermal synthesis approach would be one of the possible routes for preparing the two new compounds, as demonstrated for already known ones such as tavorite LiMgSO 4 F 37 , and structure-isotopic compounds such as LiFeSO 4 F 36 , LiFePO 4 F 50 , and LiTiPO 4 F 38,50 .
From above results, we have shown that the present DFT-coupled Bayesian optimization approach with additive structure can be applied for quinary systems and when an initial crystal structure type is provided. However, the need for an input structure means that novel compounds with new crystal structures are unsearchable. Nevertheless, we would like to point out that there is now a rich plethora of structure prototypes that can be accessed from existing databases (for example, ICSD presently contains 9,093 structure prototypes) 51 . On another note, other state-of-the-art material search methods have been reported as well, such as crystal structure prediction (CSP) techniques based on evolutionary algorithm 52 . CSP approaches do not require an input structure (the initial atomic arrangement is usually set randomly) but they need composition and initial cell volume. These techniques are meta-heuristic and utilizes empirical rules to govern the search for ground state materials. CSP techniques need to deal with the curse of dimensionality which means that local or global minima structures becomes harder and harder to find as the number of atoms and/or element type increases 52 . Combining our approach with CSP techniques, for example for quinary systems, would be one interesting direction to pursue related to high-dimensionality material search.

Conclusion
A Bayesian-driven DFT-based screening for Li and Na ionic conductors with the tavorite structure was demonstrated using ion migration energy E ( ) b as the search criterion. The BO search method was found to be generally more efficient than random search even under a stringent condition of having a positively skewed E b sample distributions. Using the Na dataset, additive BO with a knowledge transfer setting requires only an average of ~30% search space coverage to recover the optimal compound ~90% of the time. Using the same test dataset and with a search criterion of E b < 0.3 eV, additive BO also only needs to observe ~30% of the search space on the average to find ~70% of the total desired compounds, this is twice the recovery performance for desired materials of ordinary BO and random search which can only find ~37% and ~33%, respectively. These performances are realized with the use of effective descriptors, particularly RDF features. Overall, additive modeling can be an effective approach for addressing the high-dimensionality issue in BO-based searches.

Methods
Chemical search space and crystal structure description. Tavorite-type compounds with a general formula AMXO 4 Z (A: Li, Na; M: group 2, 3, 4, 13 elements; X: group 14, 15, 16 elements; Z: F, Cl, Br, I) were targeted for the E b -based solid electrolyte screening. We note that the M-X cation pair for group 5 and group 13 elements was not included in this study. Although quinary systems have been reported with group 5 cations (e.g., with Ta 5+ and Nb 5+ in another structure type 32,53 ), group 5 and 13 pairing is highly unlikely in the tavorite structure. This unlikelihood is explained by the deviation of charge distribution for the group 5 -group 13 cation pairing case which leads to a significant destabilization of the crystal structure. The model crystal structure (P 1) is shown in Fig. 1a with the host framework comprising with MO 4 Z 2 octahedra (M 1a, 1 h; O 2i) and XO 4 tetrahedra (X 2i; O 2i). The MO 4 Z 2 octahedra are corner-linked together at their trans-Z atoms to form chains along [111]. Each oxygen atoms from these chains are in turn shared with X atoms which then assumes a tetrahedral environment. Site splitting occurs for the A atoms (2i). Overall, the search space includes LiMXO 4 F dataset taken from our previous work (13) and newly calculated datasets for LiMXO 4 (Cl/Br/I) and NaMXO 4 (F/Cl/Br/I).
DFT calculation settings. The VASP code 54 was used for DFT modeling which applies the projected augmented wave (PAW) approach 55 . The energy for exchange correlation was described in the generalized gradient approximation (GGA) with Perdew-Burke-Ernzernhof formulation for solids (PBEsol) 56 . The initial coordinate dataset for the tavorite structure was referred from available crystal information file (cif) in the Inorganic Crystal Structure Database (ICSD) 52 . With a unit cell of 16 atoms and a spin-polarized condition, a 500-eV cutoff for kinetic energy and a Monkhorst-Pack kpoint resolution of 5 × 4 × 3 were confirmed to show a total energy convergence of less than 1 meV/formula unit (fu). Database-unreported tavorite compounds were calculated using the available experimental cif data as template. The calculation for static atomic charges was based from Bader method 57 . For the dynamical charges, Born effective charge calculation was carried out 58 .
The nudged elastic band (NEB) technique was employed to calculate E b 59 . The unit cell was expanded into a 1 × 2 × 2 supercell and over-the-Brillouin-zone numerical integration was performed by Γ-point sampling. With these conditions, we point out that most of the compounds especially those in the low E b region were converged to less than 10 meV/fu (with a few compounds with E b > 1.5 eV converged to less than 30 meV/fu). After structure optimization on the initial and final state supercell models containing a single A vacancy, seven images in between for the migrating A cation were constructed by linear interpolation. The value of E b was then calculated according to the formula: where E max and E min are the maximum and minimum supercell image energies, respectively, along the migration pathway.

Material descriptor formulation, formulation of DFT-E b -based search/screening driven by BO.
Candidate material descriptors were extracted from the DFT-optimized crystal structures, their description is available in Supplementary Table S1 and Supplementary Figure 1. The resulting initial domain size of the feature space has a total of 348 descriptors. Details on the construction of additive Bayesian model are provided as well in Supplementary Information section.