High-throughput production of force-fields for solid-state electrolyte materials

Research and development for solid-state electrolyte (SSE) materials is vital due to both the range of practical usage (e.g., all solid-state rechargeable batteries, fuel cells, sensors, etc.) and fundamental fascination of fast ion transport in crystalline solids.^1,2 Ion conductivity is the most representative property for SSE materials, and numerous efforts have been devoted to discover fast ion conductors. For example, all solid-state rechargeable batteries request SSE materials with sufficiently fast ion conduction comparable to liquid electrolytes,³ and thus, sometimes candidate materials are screened in terms of their ion conductivities.^4–8 Since the candidate material space is huge, it is expected to screen the candidate materials by computational approaches that can be performed for a much lower cost in terms of both financial resources and time than by experimental approaches.

There are several computational approaches to evaluate ionic migration in solids such as molecular dynamics (MD),^9,10 activation energy calculation by using the nudged elastic band (NEB) method,^11–13 kinetic Monte Carlo (kMC),^6,14 and percolation analysis.¹⁵ Among these approaches, MD can be used without a priori assumptions about the mechanism of how and where the ions move through lattice sites or they move in a concerted way, and thus, MD has been widely used to study ionic migration in solids. MD requires an accurate description of interatomic interactions, and it determines whether or not the MD simulation is meaningful and the ionic migration mechanism and activation energy predicted from the MD simulation are reliable. Ab initio MD (AIMD) has been widely used as a tool for the study of ionic migration in solids mainly because the interatomic interactions obtained via the ab initio calculation are sufficiently accurate and reliable.^16–18 However, the application of AIMD is limited to small systems (usually less than 1000 atoms) and for a relatively short time (usually less than nano-second) owing to its high computational cost. Thus, accurate and reliable empirical force-fields (FFs) are required to perform an efficient MD simulation.

Over the past decade, a lot of machine learning (ML)-type FFs have been proposed, such as linear-regression,¹⁹ neural-network (NN),^20–23 Gaussian approximation potential (GAP),^24,25 etc. It is shown that these ML-FFs can reproduce ab initio energies and forces, and thus, we could perform the MD simulation in ab initio accuracy by using ML-FFs. However, these ML-FFs require a lot of reference data, the so-called big data, to make them learn, and the optimization and validation processes usually have to be repeated several times to avoid unphysical low energy structures that can often occur because of intrinsic weakness of extrapolation in ML-FFs. For the purpose of screening of candidate SSE materials where the speed of estimation of materials properties is of importance and rough estimation would be acceptable, making an ML-FF of ab initio accuracy in a considerable amount of time is impractical.

An approach opposite of ML-FFs is to create universal FFs that are available for a part of elements in the periodic table and combinations among them. The FF based on the idea of bond valence sum (BVS) proposed by Adams and Rao^6,19,26 is a successful example of this kind. The BVS-FF can be applied to many compounds consisting of oxygen and cationic ions against oxygen,¹⁹ and it is extended to other anionic species.^6,26 Although the BVS-FF is successful on many materials, it still needs fine tuning of some parameters such as bond softness parameter when applying to a certain material,²⁶ and thus it is not obvious whether the BVS-FF can be available to specific candidate materials of interest. Thus, there is still a demand for creating empirical FFs of SSE materials with sufficient accuracy for the purpose of quickly evaluating key properties such as ion conductivity.

There is also another approach, called Iterative Boltzmann Inversion (IBI),²⁷ that uses a radial distribution function (RDF) as target data for the potential between coarse-grained particles and unit atoms in the field of polymer simulation. In the IBI, a pair potential at a distance r is iteratively updated so that the produced RDF at distance r matches the reference RDF at the same r. Since it uses a flexible tabulated form for the potential, the produced potential can well reproduce the reproduced RDF, which is the property of interest in the polymer simulation. Despite its success in the field of polymer simulation, to the best of our knowledge, the IBI has not been applied to interatomic FFs.

In this paper, we propose an automatic and high-throughput procedure of producing empirical FFs for SSE materials using radial and angular distribution functions (RDF and ADF) and equilibrium volume obtained via the AIMD simulation as reference data, partly inspired by the IBI. The FF form is similar to that of BVS-FF except that it includes a three-body angular-dependent term between selected bonds. FF parameters are optimized by using a nature-inspired metaheuristic algorithm with adaptive update of parameter search space. The detailed optimization procedure is described in Sec. II, and as a demonstration, the FF parameters of two SSE materials and their properties are shown in Sec. III.

We employ an empirical form of FF similar to that of BVS-FF^6,19 for pair interactions, that is, the screened Coulomb potential for the pairs having charges of the same signs (cation–cation and anion–anion) and Morse potential for pairs with charges of opposite signs (cation–anion). The screened Coulomb interaction is written using the complementary error function erfc(·) as

where q_A and q_B are charges of ions A and B, respectively. ρ_AB is the screening length, which is defined as a sum of effective radii of ions, ρ_AB = r_A + r_B. ε₀ is the permittivity of vacuum. Covalent-like interactions between the nearby cation and anion are described using the Morse potential,

where D_AB, α_AB, and s_AB are optimizing parameters between the interacting pair A–B. Both the Coulomb and Morse potentials are truncated at a certain cutoff radius $rc(2b)$⁠, and the following smoothing modification is applied to avoid discontinuities of potential energy and forces at $rc(2b)$⁠,

In addition to the above pair potentials, an angular potential is also applied to angles between certain selected bonds to aid forming framework structures such as tetrahedra or octahedra that are often seen in SSE materials. We adopt the Stillinger–Weber (SW)-type angular (three-body) potential,²⁸ 

where r₁ and r₂ are bond distances between the A–B and B–C species, respectively, and θ is the angle between the 1 and 2 bonds. $rc(3b)$ is the cutoff radius for this angular potential, which is usually shorter than that for two-body potentials in order that it takes only first nearest-neighbor atoms into account in the angular potential. In this paper, $rc(2b)$ and $rc(3b)$ are fixed to 6.0 Å and 3.5 Å, respectively.

There are other choices of FF function form for SSE materials, such as the full Coulomb potential including the attraction and dispersion potential. In this paper, we employ the above FF form and are not including attractive Coulomb and dispersion terms because this FF form is shown to be applicable to many SSE materials in the former studies,¹⁹ and the focus of this study is the speed of production of sufficiently accurate FFs, which usually demands the FF form as simple as possible. Results in Sec. III and some tests (not shown) indicate that this simple FF form works well for SSE materials. Since the optimization procedure described in Sec. II B and II C is applicable to any function form, incorporation of additional terms into the FF is straightforward.

In this paper, we employ RDFs, ADFs, and equilibrium volume as the reference data for optimization. The loss function for optimization is defined as a sum of squared differences of RDF g(r), ADF h(θ), and equilibrium volumes V obtained from FFMD and AIMD as

where p and t stand for the pair and triple of species, respectively, Np and Nt being the numbers of pairs and triples that are taken into account, and NR and NA being the numbers of sampling points in radial and angular distribution functions. In this study, N_R = 100 and N_A = 180 are employed.

Since the loss function, Eq. (5), is defined with RDF, ADF, and equilibrium volume that are not analytically differentiable by FF parameters, we adopt a nature-inspired metaheuristic optimization algorithm called cuckoo search (CS)^29,30 that does not require derivatives of $L$ with respect to FF parameters. The detailed description of the CS algorithm can be found in the literature.^29,30 The advantage of the CS algorithm is that it usually outperforms other metaheuristic methods such as the genetic algorithm (GA) and the particle swarm optimization (PSO) with only two hyper-parameters: the number of nests in a generation N_n and the fraction f_a of the worst nests that are to be abandoned in the next generation.²⁹ We fixed these hyper-parameters as N_n = 20 and f_a = 0.25 in all the cases performed in this study.

Even if the optimization algorithm performs well, it becomes difficult to find a good solution when the search space is too large. Thus, it is common to set lower and upper bounds to the optimizing parameters beforehand, update the bounds manually by looking at the result of optimization, and repeat this optimization and update until it reaches convergence criteria. Since this manual updating process slows down the throughput of FF production, we automate the update process of the bounds of parameters according to the following algorithm:

1. Initialize lower and upper bounds of optimizing parameters manually at the beginning of optimization.

2. Ranking: In every N steps, rank all the candidate parameter sets ever appeared with respect to $Lk$⁠, where k is a parameter-set ID, and take the top M parameter sets.

3. Shrinking: Extract the minimum and maximum values of each parameter from the M parameter sets and update the lower and upper bounds of each parameter by replacing them with the extracted minimum and the maximum values.

4. Centering the current best: Expand the lower or upper bounds so that the parameter value of the current best set is located at the center of the lower and upper bounds.

In this study, we chose N = 10 and M = 100 after some trials, and these values work well for all the cases shown in this paper. All the computations including FF calculations and CS optimization except AIMD were performed using our open-source program, nap.³¹ 

The AIMD calculation based on the density functional theory (DFT) is employed to obtain reference RDFs, ADFs, and equilibrium volume at a certain temperature. All the AIMD were performed using the Vienna Ab initio Simulation Package (VASP),³² which is based on the projected augmented wave method (PAW).^33,34 The exchange–correlation functionals are described by the generalized gradient approximation (GGA) with the Perdew–Burke–Ernzerhof format modified for solid materials (PBEsol).³⁵ The reference RDFs, ADFs, and equilibrium volume were obtained by averaging those of 20 configurations extracted from NσT-ensemble AIMD snapshots after equilibration at a certain temperature. In order to avoid the spiky shape of distribution functions, both RDF and ADF are smeared with the Gaussian function with the widths σ_RDF = 0.1 Å and σ_ADF = 2°, respectively. In most cases, the equilibration and sampling were performed with 1000 steps of AIMD, except the case that the AIMD result of our previous study was used. The equilibrium volume can converge under 1% deviation within 1000 steps. In addition, RDF and ADF seem to converge by averaging over about 50 atoms (see Sec. I of the supplementary material for more detail), and thus, we can obtain well converged reference RDF and ADF by averaging over 20 snapshots from AIMD including enough atoms per species. The same system (atomic configuration and initial cell vectors) and the same procedure to generate RDF, ADF, and volume are used in the parameter-optimization process.

Na super ionic conductor (NASICON)-type structure^36,37 is an important class of candidate SSE materials, and thus, a number of efforts have been taken to elucidate the mechanism of fast ion conduction in the structure. Among NASICON-type oxides, $LiZr2(PO4)3$ (LZP) shows relatively high lithium conductivity and has recently attracted considerable attention due to reported stable electrochemical cycles for all solid-state rechargeable Li metal battery.³⁸ The report and our previous studies^17,39 indicated the redox stability against the Li metal for LZP derivatives as well as relatively high ionic conduction of Li. However, since the ionic conductivity of Li (∼10⁻⁵ S cm⁻¹ at room temperature for Ca-doped compounds) is smaller by an order of magnitude than the other NASICON-type oxides such as LiTi₂(PO₄) and its relatives being reactive with the Li metal,^40,41 or garnet-type Li₇La₃Zr₂O₁₂-related compounds known as inert to Li metal contact,⁴² further improvement of Li ion conductivities is required.

In order to construct an FF for LZP, we apply the proposed FF production procedure and optimize its parameters to RDFs, ADFs, and equilibrium volume obtained from AIMD at 1073 K performed in our previous work.³⁹ In the AIMD, the system with 288 atoms and cell vectors of lengths (12.7 Å, 18.0 Å, and 17.9 Å) were employed. In addition, the energy cutoff for the expansion of wave functions was 500 eV, and 1 × 1 × 1 k-point (Γ-point) was used. Figure 1(a) shows the loss function values, $L$⁠, of all the parameter sets appeared during the optimization as a function of CS generation. The plot shows that the CS reaches a sufficiently small loss function value (⁠$L=0.015$⁠) within 100 generations with balancing intensification and diversification. Contributions of $LR$⁠, $LA$⁠, and $LV$ at optimum are 0.0113, 0.0037, and less than 0.000 05, respectively, indicating that the optimization effort is mainly devoted to reduce $LR$ and $LA$⁠. Figure 1(b) shows all the parameter sets projected onto two-dimensional space spanned by D_LiO and D_ZrO with points colored by $log10(L)$⁠. It is clear that the best parameter set (D_LiO, D_ZrO) = (3.273, 4.513) is outside the initial search space (parameter ranges) shown as a red broken square in Fig. 1(b), and after several updates of the search space during the optimization, the final search space covers the best parameter set. In addition, the density plots along D_LiO and D_ZrO show that there are two peaks; the smaller one in the initial search space and the larger one in the final search space. This indicates that the present method of automatic update of search space succeeds to get out from a local minimum in the initial search space and finds better parameter sets outside the initial search space, and at the same time, it intensely searches the best point in the search space.

The RDFs and ADFs obtained using the optimum parameter set are shown in Figs. 1(c) and 1(d), respectively. The FF obtained by the present procedure well reproduces RDFs and ADFs of almost all the pairs and angles, except that the difference of Li–Li RDF between FFMD and AIMD is not negligible around 3 Å. This may be caused by the fact that there are several crystallographically identical sites for the Li-ion, called 36f Wyckoff positions,⁴³ most of which are vacant, and there are a lot of Li configurations in which a part of the 36f sites are occupied. This difference in RDF would not be a significant deficient of the current FF as shown below. Comparison of bulk modulus between the DFT and present FF calculations was made using Murnaghan’s equation of state⁴⁴ to evaluate the validity of optimized FF parameters. The bulk moduli estimated using DFT and the present FF are 73.4 GPa and 87.4 GPa, respectively, showing reasonable accordance (see Sec. V of the supplementary material for more details of bulk modulus calculation). Since the bulk modulus is not referred to in the optimization of FF parameters, the above accordance supports well the validity of FF parameters.

We performed NVT-MD simulation runs of the LZP system using the present FF with parameters listed in Table I at several temperatures from 500 K to 1100 K and observed diffusion coefficients of Li-ion D_Li(T) from mean squared displacements (MSDs) of Li ions. Analyzing Li-ion probability density during the MD simulation and average distance of Li ions from the closest 6b site, it can be concluded that Li ions migrate the three-dimensional pathway along the line between the nearest-neighbor 6b sites through the 18e site, which is consistent with AIMD¹⁷ (see Sec. VI of the supplementary material for more detail). Figure 2(a) shows the MSDs of Li, Zr, P, and O ions in the MD simulation run at 1000 K, and it can be confirmed that only Li ions migrate and other ions do not, which is in agreement with the AIMD simulation.³⁹ From the Arrhenius plot shown in Fig. 2(b), we obtain an activation energy of Li-ion diffusion of 0.39 eV, which agrees well with an experimental value, 0.39 eV,⁴⁵ and an AIMD result, 0.43 eV.³⁹ 

The training data used to optimize parameters are RDFs, ADFs, and equilibrium volume obtained from a relatively short period of AIMD at a certain temperature in which atoms vibrate around their equilibrium positions. Although the training data may not include information of saddle-point structures of ion migrations, the FFs produced by the present method well reproduce the activation energy of migrating ion and stable framework structures made of non-migrating ions. This is presumably because the present FF model, consisting of the screened Coulomb, Morse, and angular potentials, is appropriate for describing the interactions between ions in SSE materials, and it holds for saddle-point structures of ion migrations. We have investigated how much the FF properties depend on the temperature at which the AIMD is performed and concluded that the diffusivity of migrating ions is not very sensitive to the temperature of the reference AIMD (see Sec. IV of the supplementary material). The angular potential may not be necessary, and FFs made of only two-body potentials could achieve a similar level of accuracy that the present FFs have. However, since the interactions between cations and anions have some covalent-bond character that has a tendency to make certain polyhedra such as PO₄, it is natural to introduce the angular terms between certain pairs of bonds. In addition, the matching of ADF between DFT and FF is improved by adding the angular terms, as shown in Fig. 1(d) and Fig. S2 of the supplementary material. Adding the angular terms to all the combinations of bonds would be redundant, and some angular terms may deteriorate the reproducibility of ion conductivity. Thus, we employ angular terms for angles between the cation–anion bonds centered at cations supposed not to migrate; for example, O–P–O and O–Zr–O are employed, but O–Li–O and P–O–P are not, in this case (see Sec. III of the supplementary material for more detailed analyses).

It is worth investigating the transferability of FFs to other systems of different structures and/or compositions. We adopt a lithium-stuffed garnet known as Li₇La₃Zr₂O₁₂ (LLZ), which is also one of the promising candidates for the SSE material. Since the peak positions of Zr–O RDFs and the shape of O–Zr–O ADF between LZP and LLZ are very similar, the FF parameters optimized to the LZP system are expected to be applicable to the garnet LLZ system. The parameters among Li, Zr, and O were fixed to those obtained for the LZP system, and only parameters related to La were optimized to RDF, ADF, and equilibrium volume obtained from AIMD of the LLZ system. Figure 3 compares the Zr–O RDFs obtained via AIMD, FFMD employing the transferred FF, and FFMD using the parameters fully optimized to the LLZ system. The peak positions of the three RDFs agree well, but the height and width of the first peak of the transferred FF differ from those of DFT and fully optimized FFMD. This indicates that the transferred FF can well reproduce equilibrium volume and lattice constants, whereas it overestimates the binding between Zr–O, which might cause some deviations in other properties such as Li-ion conductivity and thermal conductivity.

Fluoride-ion battery (FIB) is one of the candidates of post Li-ion battery, in which F ions shuttle between the cathode and anode hosts. Advantage of FIB lies in extraordinary high theoretical energy density up to 5000 Wh l⁻¹ partly due to the fact that F is the most electronegative and the lightest element among halogens.^46–48 The reversible cycling of FIB was demonstrated using SSE of BaF₂-doped tysonite-type LaF₃, La_1−xBa_xF_3−x, at elevated temperature.⁴⁹ However, the F-ion conductivity of tysonite-type polycrystalline materials is <10⁻⁵ S cm⁻¹ at 373 K (an activation energy of ∼0.5 eV), which is lower than representative Li conductive ceramics, according to the recent study.⁴⁸ Hence, further improvement and precise understandings of the F-ion conductivity in tysonite-type compounds are required. Thus, fast and accurate FFMD by using DFT-derived FF parameters is an attractive approach for optimizing composition and fundamental studies on F-ion migration.

First, we performed AIMD of tysonite structure LaF₃, which contained 192 atoms in cell vectors of lengths (12.4 Å, 14.4 Å, 14.7 Å), at 900 K. The energy cutoff was set to 400 eV, and 1 × 1 × 1 k-point (Γ-point) was employed. For doped systems, dopants were introduced by replacing randomly selected four in 48 La ions, and some F ions were added/removed to maintain charge neutrality assuming that the formal charges of ions are kept. Second, parameters related to La and F were optimized using RDFs of pairs, La–La, La–F and F–F, ADF of F–La–F (the center atoms is La), and equilibrium volume as reference data. Then, the parameters related to dopant X such as pairs X–La, X–F, and X–X and angle F–X–F were optimized to reproduce RDFs and ADF related to the dopant X, while the parameters among La and F were fixed during the optimization for dopant X. Thus, the optimization is limited compared to the case of LZP where parameters are fully optimized; for example, the loss function ends up at 0.0422, and contributions of $LR$⁠, $LA$⁠, and $LV$ are 0.0186, 0.0181, and 0.0055, respectively, which are larger than those of LZP. The FF parameters obtained by the present procedure are listed in Table II.

Figures 4(a) and 4(b) show RDFs of the LaF₃ system obtained from FFMDs with different parameter sets optimized using the different definitions of valence charges for Coulomb interactions; (a) formal valence charges (q_La = 3, q_F = −1) and (b) effective charges determined by the way suggested by Adams and Rao,^6,19 where the charges of species are reduced from the formal valences depending on their principal quantum numbers; in the case of LaF₃, the effective charges become q_La = 1.612 and q_F = −0.537. As can be seen from Fig. 4(b), the first peak of F–F RDF obtained using the effective charges is much broader than that from AIMD, indicating that the repulsive interaction between F ions is too weak and presumably leading lower activation barrier of F-ion migration. On the other hand, the RDFs (both La–F and F–F) are improved when employing the formal valence charges, as shown in Fig. 4(a). Note that the first peak of F–F RDF from FFMD is still slightly broader compared to that from AIMD, indicating that the activation energy for F-ion migration will be still slightly underestimated. Figure 4(c) shows ADFs of the F–La–F angle of FFMD with the formal valence charges and AIMD, and they are in good agreement.

RDFs and ADFs for X (X = Rb, Sr, Cs, and Ba) are shown in Fig. 5 and those of FF and DFT are mostly in good agreement, although the matchings of F–F RDF and F–La–F ADF become slightly worse than the case of the pure LaF₃ system. As it is anticipated, the peak positions of RDFs for X–F have some tendencies such that group-I elements, Rb and Cs, are slightly larger than group-II elements, Sr and Ba, and group-II elements attract more F ions than group-I elements. The present FF reproduces RDFs and ADFs well and captures the tendencies properly. However, in all cases, F–F RDFs of FFMD are broader than those of AIMD, indicating the underestimation of the activation barrier for F ion migration. We note that even though the RDFs have some tendencies with respect to groups and/or periods of elements, it is hard to find any relationship between those tendencies and FF parameters listed in Table II, e.g., magnitude relationship in D_XF and α_XF is opposite between fifth and sixth periods, such as D_RbF < D_SrF but D_CsF > D_BaF and α_RbF > α_SrF but α_CsF < α_BaF.

Figure 6(a) shows dopant-density dependence of atomic volumes at 0 K for group-II dopants, Sr and Ba. Although, the FF underestimates the volume about 10% (3%–4% in lattice constants), the dependencies on dopant species and the density of the dopant are rather in good agreement between FF and DFT; for example, Ba expands the host lattice more than Sr. Dopant-density dependence of bulk moduli is shown in Fig. 5(b), and the deviations in bulk moduli are rather notable. Bulk modulus values and their dependence on density are both underestimated. However, the dependence on dopant species is correctly reproduced, e.g., as x increases, the bulk modulus decreases more for Ba than for Sr. Although these deviations could be improved by introducing the volume at 0 K and bulk modulus into the loss function in Eq. (5), as discussed in the following, the quality of the present FF is presumably sufficient for the purpose of studying the doping effects on ion conductivity and screening candidate SSE materials.

To investigate the quality of the present FF for La–X–F on F-ion migration, we performed NVT-ensemble MD simulations of non-doped tysonite LaF₃, doped La_1−xBa_xF_3−x, and La_1−xSr_xF_3−x with x = 0.083 at several temperatures. Figure 7 shows the Arrhenius plots of F-ion diffusion coefficients obtained from FFMDs as a function of 1000/T, and activation energies of F ions in LaF₃, La_1−xBa_xF_3−x, and La_1−xSr_xF_3−x are 0.61 eV, 0.41 eV, and 0.35 eV, respectively. There is a wide variety in reported activation energies of LaF₃ from 0.18 eV to 1.2 eV^48,50–52, and it is also suggested that the activation energy depends on temperature.^50,51 The activation energy of LaF₃ from 600 K to 1200 K obtained with the present FF, 0.61 eV, can be compared with 0.76 ± 0.03 eV,⁵⁰ 0.8 eV–0.85 eV⁵² in the range from 550 K to 850 K and 0.76 eV⁴⁸ at lower temperatures than 573 K, indicating that the present FF underestimates the activation energy of F-ion migration. This underestimation is probably related to the broadness of the first peak of F–F RDF as mentioned above.

The drop of the activation energy from LaF₃ to La_1−xSr_xF_3−x by doping obtained by the present FF, 0.26 eV, is comparable to that in an experimental result,⁴⁸ 0.27 eV (from 0.76 eV to 0.49 eV.) This indicates that the present FF correctly reproduces the effect of strain created by Sr ions doped in LaF₃ and its effect on the F-ion migration. The activation barrier of La_1−xBa_xF_3−x is lower than that of LaF₃ and slightly higher than that of La_1−xSr_xF_3−x, which is considered as a result of the difference of ionic radius between Ba and Sr and agrees well with the experimental result.⁴⁸ These results demonstrate that the FF for dopants whose parameters are optimized with fixed parameters between host species works well to study the doping effects on ion conductivity.

Traditionally, the FF parameters are optimized to static parameters such as lattice constants, bulk moduli, and elastic constants obtained by experiments or ab initio calculations. For complex systems involving over two elements, which is common in SSE materials, such data are sometimes unreliable or even unavailable, and thus, the force-matching method⁵³ has also been widely used for both the empirical and ML-type FFs, for example, an open-source program potfit^54,55 employs the force-matching as well as matching to energy and stress. It has been shown that the force-matching method provides accurate FFs if reference data are carefully prepared and the matching is sufficiently good. However, it is not an easy task because the appropriate reference data could be different to each structure or composition, and unlike ML-type FFs, it is common that empirical FFs cannot achieve good matching of forces. We have tested the force-matching method for several SSE materials with the same reference data as used in the present DF-matching method, i.e., forces in snapshot structures from AIMD, and in some cases, the qualities of FFs were poor, and crystal structures were not kept even at low temperature. On the other hand, the present DF-matching method performs the MD simulation during the optimization, and thus, it is guaranteed that the FFs reproduce stable crystal structures up to the temperatures used in the optimization. Since the reproduction of RDF and ADF is an indication of good prediction of interactions between atoms, the produced FFs usually achieve high accuracy not only in static properties but also in dynamic properties such as ion conductivity.

We note about the difference from the IBI.²⁷ In the IBI, a tabulated potential at distance r is iteratively updated to improve the ratio of reference RDF and the RDF obtained by the potential at the same r. On the other hand, in the present approach, the FF function form is fixed, and all the RDFs and ADFs of selected pairs and triples are optimized simultaneously because, in the case of SSE materials, interactions between species contribute to the RDFs and ADFs related to the other species in the structure. Optimizing the parameters simultaneously not only to RDFs but also to ADFs and volume is presumably making it possible to produce high-quality interatomic FFs.

As for the computational cost for the production of FFs, the most time-consuming part is to obtain reference data (RDFs, ADFs, and volume) from the AIMD simulation, which usually takes several hours to a day even by using high-performance computers. In the process of optimization of FF parameters using the CS algorithm, our current implementation is able to parallelize evaluations of nests in a generation using processor cores in a computer node. Thus, by using a workstation of over 20 cores, the optimization process takes a few hours to perform 100 generations, which is sufficient for all the cases shown in this paper. In total, a FF can be produced within a day or two using the present approach.

In conclusion, an automatic and high-throughput method to produce FFs for SSE materials is proposed. The present FF model employs screened Coulomb repulsions between cations and anions, the Morse potential between cations and anions, and the angular potential between the cation–anion bonds connected at cations that are not supposed to migrate. RDF, ADF, and equilibrium volume obtained from AIMD are used as the reference data for the optimization of FF parameters. The CS algorithm with an automatic update of search space, shrinkage, and expansion of parameter ranges works well for the purpose of automatic determination of FF parameters. FFs for LZP and La–X–F systems are created using the present method, and the activation energies of migrating ions obtained via FFMD are shown to be in good agreement with those obtained from experiments. All the FF parameters listed in the present paper were obtained through this automatic procedure, which took only a day or two including the reference-data production by AIMD. The present approach enables us to create a FF of any novel SSE material once 1000-step AIMD is performed and to perform a sufficiently accurate MD simulation of ion migration, which will greatly enhance the speed of materials exploration and screening.

The supplementary material provides the following analyses and/or descriptions: convergence of distribution functions with respect to the number of sampling atoms, ADFs of the LZP system that are not shown in the main text, sensitivity analyses on angular terms of the FF, sensitivity analysis on the reference AIMD condition, bulk modulus calculation, and the analysis of Li positions during MD in the LZP system.

The computer programs used to perform the CS algorithm and MD simulation during optimization are available in the open-source package, nap.³¹ The data required to reproduce this study such as AIMD snapshots, RDFs, ADFs and volumes of LZP and La–X–F systems are available in the supplementary material.

This work was supported in part by the “Materials Research by Information Integration” Initiative (MI2I) project of “Support Program for Starting Up Innovation Hub” from the Japan Science and Technology Agency (JST), JSPS KAKENHI (Grant Nos. JP18H01707, JP19H05815*, JP20H02436, and JP20H05290*) (*Grant-in-Aid for Scientific Research on Innovative Areas “Interface Ionics”), “Elements Strategy Initiative to Form Core Research Center” from the Ministry of Education Culture (MEXT), and funded by Toyota Motor Corporation. The authors would like to thank Dr. Hisatsugu Yamasaki and Dr. Kazuto Ide from Toyota Motor Corporation for fruitful discussion. The authors would also like to thank the Information Technology Center of Nagoya University for providing computational resources.