Strategies for fitting accurate machine-learned inter-atomic potentials for solid electrolytes

Ion transport in materials is routinely probed through several experimental techniques, which introduce variability in reported ionic diffusivities and conductivities. The computational prediction of ionic diffusivities and conductivities helps in identifying good ionic conductors, and suitable solid electrolytes (SEs), thus establishing firm structure-property relationships. Machine-learned potentials are an attractive strategy to extend the capabilities of accurate ab initio molecular dynamics (AIMD) to longer simulations for larger systems, enabling the study of ion transport at lower temperatures. However, machine-learned potentials being in their infancy, critical assessments of their predicting capabilities are rare. Here, we identified the main factors controlling the quality of a machine-learning potential based on the moment tensor potential formulation, when applied to the properties of ion transport in ionic conductors, such as SEs. Our results underline the importance of high-quality and diverse training sets required to fit moment tensor potentials. We highlight the importance of considering intrinsic defects which may occur in SEs. We demonstrate the limitations posed by short-timescale and high-temperature AIMD simulations to predict the room-temperature properties of materials.


Introduction
Rechargeable lithium (Li)-ion batteries are widely used in portable electronics and are seeing applications in powering electric vehicles. To be competitive with internal combustion * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. engines, batteries for automotive applications must display energy densities >350 Wh kg −1 and cost below $100 kWh −1 [1]. One of the major safety concerns of commercial Li-ion batteries arises from the flammable nature of the liquid organic electrolytes. Nonflammable solid electrolytes (SEs) are promising alternatives to mitigate safety issues [2,3].

Future perspectives
One of the major challenges in all-solid-state batteries is the incompatibility of the interface between the solid electrolytes (SEs) and the electrode materials, which leads to the formation of resistive or reactive decomposition products during battery cycling. These reactions drive the formation of solid-solid interphases, which may be highly resistive to ionic transport with continuous thickness growth, thus contributing to the increased impedance associated to the transport of Li-ions across the interface. Thus, much-needed work is to understand the Li-ion transport across bulk materials that may compose these complicated heterogeneous multiphased electrode-electrolyte interfaces. Simulation techniques set a firm link between phenomena occurring at the atomic scale and macroscopic observations of energy materials, i.e. SEs and electrodes. Machine-learning molecular dynamics (MD) are orders of magnitude computationally less expensive compared to ab initio MD but provide comparable accuracies. Machine-learning MD simulations based on the moment tensor potentials method represent a powerful methodology to reveal the complex nature of ion transport in SEs and their decomposition products formed upon contact with the electrode materials. assessed by several electrochemical techniques, such as galvanostatic or potentiostatic intermittent titration techniques [12,13] as well as electrochemical impedance spectroscopy [12,14]. Solid-state nuclear magnetic resonance [15], and quasi-elastic neutron scattering [16,17] are also used to measure ion diffusivities in materials. In principle, the ionic diffusivities measured by different techniques should be similar. However, the yielded values could differ by more than two orders of magnitude from different experimental techniques [12,[18][19][20]. Therefore, it is important to provide accurate calculations and reveal the ionic transport mechanisms by using computational approaches. Recently, Deng et al [21] have showcased extensive simulation of ion-transport in a popular SE, Na 1+x Zr 2 Si x P 3−x O 12 with kinetic Monte Carlo (kMC) bearing the accuracy of first-principles calculations, which can improve upon the obvious limits of ab initio molecular dynamics (AIMD) and classical molecular dynamics (MD) simulations. Albeit kMC is a powerful technique, fitting accurate lattice Hamiltonians is not straightforward.
Theoretical tools, such as AIMD, have been used extensively to describe ion transport in functional materials [22][23][24][25][26]. However, due to their high computational costs, AIMD simulations can treat at most a few hundred of atoms approaching limited simulation timescales of hundreds of picoseconds (ps), which extensively limits the required diffusion statistics associated to ion motion in materials. An alternative is classical MD simulations based on empirical force fields (FFs) which can treat much larger systems (thousands of atoms), and extend to timescales of nanoseconds (ns). FFs are usually derived from rigorous physical models and can be used across systems with similar chemical constituents, but their accuracy is also limited by their functional forms. Furthermore, the derivation of FFs is not trivial, and in some cases (e.g. charge optimized many body potential and ReaxFF) hundreds of parameters need to be determined [27][28][29].
Usually, the fitting of MTPs requires the construction of training sets and the optimization of several hyper-parameters, using which the 'learning' of the FF is achieved. Training sets can be constructed using structures from ab initio calculations, whose quality will primarily affect the accuracy of MTPs [31][32][33].
In this work, we performed AIMD simulations for generating training sets for the MTPs and explored various strategies or variables for training the MTPs. Firstly, different AIMD training sets were generated by considering several computational and physical variables. An immediate computational variable is the density functional theory (DFT) exchange and correlation (XC) functional used for the AIMD simulations, such as the standard generalized gradient approximation (GGA) (details in section 2) [37]. Using a van der Waals (vdW)-corrected XC functional, Qi et al [32] demonstrated the predicted activation energies and ionic conductivities of three 'compact' lithium superionic conductors, i.e. Li 0.33 La 0.56 TiO 3 , Li 3 YCl 6 and Li 7 P 3 S 11 are in good agreement with existing experimental values. Other computational variables are the AIMD simulation time which affects the size of the training sets and the simulation temperatures at which AIMD simulations were performed. The defect types present in the SEs and their concentrations are considered as physical variables Secondly, several MTP fitting parameters directly affect the quality of MTPs and are explained in detail in section 3.
The generated training sets were subsequently used to fit MTPs, and large-scale MD simulations were performed to predict the activation energies (E a ), ionic diffusivities, and conductivities (σ) of three topical superionic conductors, i.e. two argyrodites Li 6 PS 5 Cl and Li 6 PS 5 I and the α-phase tetragonal Na 3 PS 4 .
We find that for argyrodite Li 6 PS 5 Cl, E a predicted by MTP fitted with training data produced by a conventional GGA XC functional (PBEsol, see section 2) achieved the best agreement with existing experiments [38]. The estimated values of room temperature conductivity σ RT from MTPs fitted with different strategies are also in the range of the experimentally reported values [5,[38][39][40][41]. Our investigation demonstrated that including a sufficiently diverse set of structures bearing different local environments appears important to achieve accurate simulations.

AIMD for training sets
We used the cubic (with a F43m space group) structure of argyrodites Li 6 PS 5 X (X = Cl, I), containing 52 atoms for each structure. The initial structures were taken from Kraft et al [38]. We enumerated all symmetrically distinct orderings of the primitive cell (13 atoms). In total 26 structures were obtained for Li 6 PS 5 Cl/Li 6 PS 5 I and their total energies were evaluated using single-point DFT calculations. The lowest energy structures were fully relaxed using DFT and then used for further AIMD simulations.
In the DFT calculations, wavefunctions were described by expanding the valence electrons with plane waves together with projected augmented wave potentials (PAWs) for core electrons, as implemented in the Vienna Ab-initio Simulation Package (VASP) [42][43][44]. In the PAWs, the valence electron configurations for each element were set as follows: Li: s 1 p 0 , Na: s 1 p 0 , P: s 2 p 3 , S: s 2 p 4 , Cl: s 2 p 5 and I: s 2 p 5 . A plane wave energy cutoff of 520 eV and a k-point mesh with a k-point density of 25 Å −1 were used.
Preliminary geometry optimizations of the known structures of superionic conductors were achieved using the MITRelaxSet set implemented in pymatgen [45]. The total energy of each structure was converged to within 10 −5 eV/cell, and the geometry optimizations were assumed to converge when the change in total energy between two subsequent ionic steps was below 10 −4 eV.
After geometry optimizations, AIMD simulations to generate the training sets for the fitting of MTPs were performed using VASP. Here, the plane wave energy cutoff was reduced to 400 eV, and the total energy was integrated at the Γ-point. The canonical ensemble (i.e. at fixed number of particles, volume, and temperatures, NVT) was achieved using Nosé-Hoover thermostat and a timestep of 2 fs [46,47].
We tested several XC functionals while performing AIMD simulations. The 'standard' GGA as parameterized by Perdew, Burke, and Ernzerhof (PBE) [37] and its revision for solids PBEsol [48]. optB88 van der Waals (vdW) [49] DFT functional was used in a recent study [32] to calculate ion transport properties in Li-based superionic conductors. The strongly constrained and appropriately normed (SCAN) [50] is a recent meta-generalized-gradient approximation and was shown to be accurate in the prediction of geometries, formation energies, and migration barriers of diversely bonded molecules and materials (including covalent, metallic, ionic, hydrogen and vdW bonds) [51][52][53] SCAN + rVV10 [54] is a versatile vdW density functional that combines the SCAN meta-GGA XC functional with the rVV10 non-local vdW correlation functional. R2SCAN [55] is a regularized-restored SCAN functional that significantly improves the numerical efficiency and accuracy under low-cost computational settings as compared to its parent SCAN functional.
For argyrodites Li 6 PS 5 X (X = Cl, I), we performed AIMD simulations at 500 K, 1000 K, 1500 K, and 2000 K, each for 8 ps (preceded by a temperature ramping of 2 ps), which in total resulted in training sets containing 16 000 snapshots. To test the strategy 'Simulation time' for argyrodites, we performed AIMD for 4 ps at each temperature and included 8000 snapshots in the training set. While for α-Na 3 PS 4 , the temperatures chosen ranged from 400 K to 1600 K with a step of 400 K. The default training sets for Na 3 PS 4 contain 16 000 configurations, 4000 for each temperature. To test the strategy 'Simulation time' for Na 3 PS 4 , we performed AIMD for 40 ps at each temperature and included 80 000 snapshots in the training set.
To study the defect-mediated diffusion events in Li 6 PS 5 Cl, four types of defects were considered, i.e. Li interstitials Li i , which are the intrinsic defects in argyrodites, antisites Cl S and S Cl , and Li vacancies V Li . For antisites Cl S , from the conventional cell, we replaced two S atoms with two Cl atoms (2% defect concentration), and two Li atoms were removed from the cubic unit cell to maintain charge neutrality. In the case of the S Cl antisite defect, starting from the conventional cell, we replaced two Cl atoms with two S atoms (50% defect concentration), and two Li atoms were added in the cubic unit cell to maintain charge neutrality. In the case of V Li , a 2 × 1 × 1 supercell containing one Li + vacancy was created, corresponding to a vacancy concentration of ∼2%. Since limited diffusion events were found in the pristine tetragonal α-Na 3 PS 4 , we created a 2 × 2 × 2 supercell model containing one Na + vacancy, corresponding to a concentration of 2%. For both, Li 6 PS 5 X (X = Cl, I) and Na 3 PS 4 , Li + /Na + vacancies were introduced by removing Li/Na atoms and compensating with a uniform (jellium) charge background.

Moment tensor potential-molecular dynamics
MTPs of the superionic conductors investigated in this work were developed using the machine learning of the inter-atomic potentials package [33]. The MTP-MD simulations were performed using LAMMPS [56]. Nosé-Hoover thermostat was used to simulate the canonical ensemble (NVT) [46,47]. Long MD simulations were carried out for at least 10 ns with a short timestep of 1 fs, preceded by a temperature ramping for 10 ps (to reach each target temperature) and an equilibration period of 1 ns. To investigate the effect of the concentration of defects on the ion diffusion of Li 6 PS 5 Cl, we created supercells with different sizes. Table 1 details the structure models investigated. In the case of Cl S Concentration-1, we replaced two S atoms with two Cl atoms, and two Li atoms were removed to keep charge neutral in the conventional cell (1 × 1 × 1). For other cases, we replaced one S atom with one Cl atom, and one Li atom was removed to keep the charge neutral in the respective cell sizes. For S Cl , we only considered a defect concentration of 1.64 × 10 19 cm −3 , where we replaced one Cl atom with one S atom and added a Li atom to keep charge neutral in a 4 × 4 × 4 supercell (table 1). For V Li+ , we removed one Li atom as the size of the supercell model was increased.
where r i (t) is the displacement of the ith ion at time t, N is the number of diffusing ions, t is time and d is the dimensionality of the system. E a is the Li/Na ion migration energy and can be derived from the Arrhenius relation of equation (2), D 0 is the ion diffusivity at infinite temperature and k B is the Boltzmann constant. The ionic conductivity at temperature T is obtained from the Nernst-Einstein equation (3).
where ρ is the molar density of diffusing ions in the unit cell, z = 1 is the charge of Li and Na ions, and F and R are the Faraday constant and gas constant, respectively.

The distinct van Hove correlation function.
The distinct part of the van Hove correlation functions (G d ) of superionic conductors was calculated by using the trajectories from MTP-MD simulations as shown in equation (4).
where r i (0) is the position of the ith ion at initial time t = 0, r j (t) is the position of the jth ion at time t, N d is the number of diffusing Li/Na ions in the unit cell, and r is the radial distance between ions. The average number density ρ is used to normalize G d so that G d → 1 when r → 1.

Generalized phonon density of states.
To calculate the partial atomistic vibrational density of states, a Fourier transform of the velocity autocorrelation function [57] of the trajectories of individual atoms (Li, Na, P, S, Cl, I) was performed. The partial vibrational density of states was then neutron-weighted to obtain the generalized phonon density of states (GDOS). GDOS involves a weighting of the scattering species with scattering powers σ/M, where σ is the neutron scattering cross section in barn (1 barn ∼10 −24 cm 2 ) and M is mass of an ion in the atomic mass unit or amu (1 amu ∼0.66 × 10 −24 g). The scattering powers in barn·amu −1 of Li, Na, P, S, and Cl, I used were 0.197, 0.143, 0.107, 0.032, 0.474, and 0.030, respectively. Figure 1 illustrates the strategies and the variables associated with fitting MTPs for 'large' scale MD simulations.

Strategies for constructing moment tensor potentials
In MTP, the training sets are usually derived from AIMD simulations, and their qualities largely depend on several computational and physical variables. Some important computational variables are (a) the choice of DFT XC functional, (b) the simulation temperatures (controlling the probabilities of migration events) at which the AIMDs are executed, and (c) the simulation time of each AIMD, relating to the size of training sets. Different types of defects existing in the ionic conductors discussed above have been considered important physical variables. The computational variables are explained in greater detail in the Methods section.
Once the preliminary AIMD simulations are performed (step no. 1 in figure 1) the fitting of the MTP potential starts. MTPs are fitted following the force-matching minimization strategy of Ercolessi and Adams [58], which minimizes the residual between the forces predicted by AIMD simulations and those obtained from a trained MTP. At the stage of MTP training, several parameters need to be chosen carefully. The lev max is used to determine the completeness of the basis set defining the functional form of MTP [33]. Using too large values of lev max can easily lead to overfitting. For example, while training an MTP for argyrodite and Na 3 PS 4 , a lev max of 12 provides adequate levels of fitting and validation errors in energy (<15 meV atom −1 ) and forces (<140 meV Å −1 ), as summarized in tables S1-S4 of the supplementary information (SI). Another important parameter required in the fitting of MTPs is the radius cutoff (R cut ) used to ensure the smooth behavior of MTPs when atoms leave or enter the interaction neighborhood [33]. Here, we chose R cut to be 5 Å, which has been used in previous reports [32,33,[59][60][61]. While fitting the MTPs, a reasonable choice of weights associated to the energy and force terms is also needed. In this work, we used a ratio of weights of 10:1 for energies to forces, similar to our previous work [61]. The original training set containing AIMD snapshots is randomly split into two parts-a 'training set' and a 'testing set' with a ratio of 90:10. The MTPs are then trained on the 'training sets' and the 'testing sets' are used to validate the MTPs.
The fitted MTPs were then used to perform large-scale MD simulations. To evaluate the performance of the MTPs generated using different strategies, we compare the predicted tracer diffusivities (D * ), ion conductivities (σ), and activation energies (E a ) of superionic conductors against experimentally reported values.

Evaluation of MTPs
The crystal structures of the superionic conductors considered in this work are shown in figure S1 of the SI. In figure S2 and table S5, we show the predicted lattice parameters and volumes of the superionic conductors by using different DFT XC functionals. For argyrodite Li 6 PS 5 Cl, the volume predicted by optB88 functional shows the smallest discrepancy (∼0.04%) concerning the experimental data [38,62]. In contrast, the use of meta-GGA XC functionals, in this case, SCAN appears to overestimate the volume of Li 6 PS 5 Cl (∼10.3%), and the predicted error decreased to ∼8.4% with the addition of vdW interactions. Unsurprisingly, the performance of R2SCAN (∼10.2%) appears similar to that of SCAN. In the case of Li 6 PS 5 I, optB88 achieved the best agreement with the experimental value, with a slight overestimation of ∼0.15%.
We first look at how the training sets generated using different DFT XC functionals affect the predictive accuracy of MTPs in evaluating the ion transport properties in the superionic conductors. For clarity, we will label the MTPs fitted using different DFT XC functionals as MTP@PBE, MTP@PBEsol, MTP@optB88, etc. Figure 2 shows the variability of the activation energies E a obtained via MTPs trained with different computational variables (see figure 1). In figure 2(a), E a of Li 6 [38].
The experimentally reported ionic conductivities of Li 6 PS 5 Cl at room temperature range from 0.01 mS cm −1 to a few mS cm −1 [5,[38][39][40][41]. From our MTP-MD simulations (table S6 of SI), σ RT predicted by MTP@PBEsol was the lowest (∼0.070 mS cm −1 ), while σ RT predicted by other MTPs (trained with different DFT XC functional) assumed higher values (∼0.145-0.24 mS cm −1 ). E a and σ RT predicted from MTPs fitted using different DFT XC functionals for Li 6 PS 5 I are shown in table S7 of SI. The predicted values of E a are higher than the experimental value [38] by more than ∼150 meV ( figure S3(a)). The estimated values of σ RT are in the order of 10 −5 -10 −4 mS cm −1 , which is nearly one order of magnitude lower than the experimentally reported values (table S7) [38].
In the case of the pristine (without native defects, e.g. Na vacancies) tetragonal α-Na 3 PS 4 , we found limited diffusion events due to which it is hard to derive accurate values of diffusivity from the computed mean square displacement (MSD) plots. Thus we studied the Na + diffusion in Na 3 PS 4 , containing ∼2% of Na + vacancies. The predicted values of E a by MTPs fitted with different functionals are quite similar (see figure S4 zoom-in plot) but are lower than the experimental data by more than 0.2 eV ( figure S4). This overestimation of Na-ion diffusivities is amplified even further in the Na-ion conductivities (σ), which are three orders of magnitude higher than the experimental value (table S8) [63].
Besides the choice of DFT XC functional, it is also important to investigate the effects of other computational variables, such as the simulation temperatures and timescales used to generate the AIMD training sets. To evaluate these variables, we keep the DFT XC functional as PBEsol. We discuss the effect of simulation temperatures on MTP predictions. Our default training sets include snapshots generated from AIMD simulations performed at elevated temperatures, i.e.  in panel (a)). 500 K, 1000 K, 1500 K, and 2000 K for argyrodites; and 400 K, 800 K, 1200 K, and 1600 K for Na 3 PS 4 . In the case of argyrodites, we removed the low-temperature snapshots (500 K) from the training sets, and the fitted MTPs are labeled as high-T-MTP. For both Li 6 PS 5 Cl and Li 6 PS 5 I, the predicted values of E a using high-T-MTP training sets are lower than the results from MTP@PBEsol by ∼70 meV (figure 2(b) and table S6 in SI) and ∼30 meV (figure S3(b) and table S7 in SI) respectively. A decrease in values of E a may also lead to an overestimation of the ionic conductivities σ RT (tables S6 and S7). For Na 3 PS 4 , since the predicted E a were lower than the experimental values, regardless of the XC functional used, we removed the high-temperature snapshots (1600 K) from the training sets, to examine if the discrepancy between the predicted and experimentally reported values of E a could be reconciled. Unfortunately, the predicted E a and σ RT appear similar to the values obtained using MTP@PBEsol (table S8).
We then consider the effect of the timescale of AIMD simulations on the accuracy of MTPs predictions. Longer AIMD simulations include additional trajectories in the training sets, which in turn lead to increased computational costs, but also a richer dataset for the training procedure. For both argyrodites and Na 3 PS 4 , the default training sets include 16 000 snapshots (4000 snapshots for each temperature, see Methods section). For argyrodites, given that the E a predicted by MTP fitted using the default training sets show reasonable agreement with experimental values, we investigated whether the size of the training sets can be further reduced, while still maintaining the predictive capabilities of the MTP. We used training sets containing 8000 snapshots (2000 for each temperature). The predicted values of E a appear remarkably different from the experimental values (figure 2(b) and table S6 for Li 6 PS 5 Cl, figure S3 and table S7 for Li 6 PS 5 I), indicating that a reduction in the size of the training sets cannot guarantee similar accuracy in the prediction of E a . In the case of Na 3 PS 4 , we performed appreciably longer AIMD simulations and included more snapshots in the training sets (∼80 000 in total, see Methods section). However, using enlarged training sets, no significant change was observed in the predicted E a , with the new predictions falling within the error bar of MTP@PBEsol (see figure S4).
The calculated tracer diffusivities D * , and conductivities σ of Li + in argyrodite Li 6 PS 5 Cl are shown in figure 3. The values of E a and σ RT are summarized in table S6.
Values of E a are derived from the Arrhenius plots, whereas room temperature conductivities (σ RT ) are extrapolated from high-temperature data. At low temperatures the number of rare events leading to Li + migration is scarce, and the small time propagation (hundreds of ns) of our MTPs is not sufficient to extract (or extrapolate) reliable low-temperature conductivities and diffusivities [21]. At low temperatures, the number of rare events leading to Li + migration may be limited, which makes the direct calculation of reliable low-temperature diffusivities and conductivities difficult.

Effect of native defects on ionic diffusivity
Native defects exist in all the materials with sizeable effects on the intrinsic electronic and the ionic conductivities of SEs [64][65][66]. However, experimentally it remains a major challenge to determine the defect types and their concentrations. A recent study by Gorai et al [64] investigated the chemistry of point defects (vacancies, antisites, and interstitials) in argyrodites Li 6 PS 5 X (X = Cl and I) and Na 3 PS 4 , by employing a combination of DFT and semi-empirical models. The derived electronic conductivities of the superionic conductors achieved quantitative agreements with the experiments [64].
Here, we further investigate the effects of these point defects on the ion transport of argyrodite Li 6 PS 5 Cl. Since the reported defect concentrations are dilute [64], sufficiently large supercell models are required, which may be inaccessible via AIMD simulations. With MTPs, we can perform largescale simulations including thousands of atoms, thus capturing scenarios of the dilute defect concentrations. Note, the defect concentrations used in AIMD simulations are much higher than those used in the MTP-MD. The details of the supercell models used in both AIMD and MTP-MD can be found in the Methods section. Argyrodite Li 6 PS 5 Cl shows intrinsic Li interstitial defects, and we analyzed their activation energies in figure 2. To distinguish from other types of defects, we use a simplified Kröger-Vink notation with Li interstitials indicated as Li i , anion antisites as S Cl and Cl S , and Li + vacancies as V Li [64]. The MTPs fitted using training sets generated from AIMD simulations for structures containing specific types of defects, are labeled as MTP Li i , MTP S Cl , MTP ClS , and MTP V Li , respectively. In these simulations, the DFT XC functional used is PBEsol.
The calculated Li + diffusivities and conductivities of Li 6 PS 5 Cl containing different types of point defects in the limit of dilute concentration are shown in figure S6. The activation energies derived from the Arrhenius plots are shown in table S9. We find that the predicted E a of Li i , Cl S and V Li are similar and are close to the experimental value reported by Kraft et al [38], suggesting the coexistence of these defects in the material together with the intrinsic interstitial defects. The estimated σ RT of these defects are also in the same order of magnitude (table S9). To investigate whether the MTP trained for Li 6 PS 5 Cl with only intrinsic defects Li i (MTP Li i ) can be extended to other defect types, we performed large-scale MD simulations for Li 6 PS 5 Cl models containing S Cl , Cl S and V Li , as shown in figure 4.
In the case of Cl S and V Li , the predicted E a using MTP Li i is similar to that from MTP ClS and MTP V Li , respectively. In contrast, by introducing S Cl in Li 6 PS 5 Cl, E a predicted by MTP Li i and MTP S Cl show appreciable discrepancy of ∼0.1 eV.
To further determine the effect of the defect concentrations, we performed MTP-MD for V Li and Cl S at different defect concentrations (table 1) using both MTP Li i , and MTPs fitted using training sets containing V Li and Cl S , respectively. The calculated Li + diffusivities and conductivities are shown in figures S7-S10 and tables S10-S11 of the SI. The derived E a as a function of defect concentrations is shown in figure S5. In the case of V Li , the dilute regime spans from 1.31 × 10 20 cm −3 to 1.05 × 10 18 cm −3 , consistent with the fact that Li + vacancies are among the lowest formation energy defects in Li 6

Lattice dynamics of Li 6 PS 5 Cl and correlation effects
It is of great interest to understand the origin of ion migration in ionic conductors. Experimentally, the ion dynamics in crystals can be revealed from techniques, such as inelastic neutron scattering spectroscopy (INS), and complementary to the more  [64]. Exp. is for the experimental data, and others in the legend refers to other defects. common Raman spectroscopy. It is established that INS is a more sensitive technique than Raman spectroscopy for detecting lower energy modes [62,67]. INS measurements can be reproduced and predicted by appropriately weighted phonon generalized density of states (GDOS) [62]. The accessibility of long MD trajectories via MTP-based simulations enables us to predict the GDOS, and thus establish a robust link between structure-property relationship [62].
Here, we calculated the GDOS of argyrodites Li 6 PS 5 X (X = Cl, I) and α-Na 3 PS 4 at different temperatures by taking the Fourier transform of the velocity autocorrelation functions [57], as obtained from our MTP PBEsol -MD simulations (see Methods). GDOSs extracted from MD simulations contain anharmonic effects, which often involve vibrational modes of species with different masses. However, this analysis is unable to probe into the type of vibrational modes, which can be accessed through DFT simulations or by direct fingerprinting of Raman spectra.
The calculated GDOS of Li 6 PS 5 Cl at 300 K and 500 K are shown in figure 5. The Li partial GDOSs contribute to the entire energy range (0-200 meV), while Cl atoms show contributions to the peak at ∼15 meV. From the zoom-in plots, the P and S partial GDOS show similar features (figure S13 of SI) motion. We observe that a temperature increase promotes the broadening of GDOS and an increase in the height of the low-energy peaks. This indicates an enhancement of the Li + motions at higher temperatures. The peak positions did not show a significant variation at elevated temperatures, which is an indication of the intrinsic softening of Li 6 PS 5 Cl [68]. Our results are both qualitatively and quantitatively consistent with recent INS measurements and MD simulations on other thiophosphate systems [68].
The calculated GDOS of Li 6 PS 5 I at 300 K and 500 K are shown in figures S14 and S15, exhibiting similar features to the GDOS of Li 6 PS 5 Cl. The GDOS of α-Na 3 PS 4 from our MTP-MD simulations (figure S16) also agree with previous INS measurements and AIMD simulations [62].
The van Hove correlation function is a useful analysis for investigating the correlated motions of ions in crystalline ion conductors [69,70], which can be divided into the self part (G S ) and distinct part (G d ) [69,70]. In particular, G d is the probability density of finding at time t a different moving ion at a distance r from a reference ion that was initially located at the origin at time t = 0. Thus, G d can access the collective dynamics of the moving ions. The calculated G d for Li + in argyrodites Li 6 PS 5 X (X = Cl, I) and Na + in Na 3 PS 4 are shown in figures 6 and S17-18, respectively.
In figures 6(a) and (b), we find large amplitudes of G d for r < 2 Å, which is a strong indication of collective motions in both Li 6 PS 5 Cl and Li 6 PS 5 I. By inspecting the initial 10 ps of the MD simulations (see figure S17 ), a larger peak of G d was identified in Li 6 PS 5 Cl, but not in Li 6 PS 5 I. This suggests that the extent of collective motions of Li + are stronger in Li 6 PS 5 Cl than in Li 6 PS 5 I, which explains faster Li + diffusion in Li 6 PS 5 Cl. We also find pronounced signatures of G d at the proximity of r = 0 Å in Na 3 PS 4 (figure S18), suggesting that the motions of Na + are highly correlated in this material [62].

Discussion
Ion transport in materials is routinely determined through several experimental techniques, which introduce a source of variability in the values of ionic diffusivities and conductivities. Likewise, a certain degree of variability in these values could also arise from simulations with different methodologies used for prediction. Here, we analyzed and revealed the main variables involved in the construction of reliable MLPs, in the flavor of MTPs. We focused on two topical SEs, which are the argyrodites Li 6 PS 5 X (with X = Cl and I) and Na 3 PS 4 .
Specifically, we identified, three main classes of variables affecting the predictive power of MTPs (see figure 1): (a) physical variables, (b) computational variables, and (c) MTP fitting parameters. The physical variables include the types of defects and their concentrations in SEs. The accuracy of MTPs depends on several computational variables, such as the types of DFT XC functional, the simulation time, and temperature, respectively, in preparing the AIMD training sets. A number of MTP-specific parameters also require attention during the fitting of reliable MTPs. We have commented on the importance of identifying appropriate values of the hyper-parameter lev max . The values of lev max that are too small may lead to inaccurate results, whereas larger values will cause overfitting. From our experience, in complex mixed polyanion superionic conductors, such as NaSICON (Na 1+x Zr 2 Si x P 3−x O 12 ) an optimization of lev max is necessary to fit accurate potentials.
The sample preparations, such as the synthesis conditions, show up as variations in the types of defects and their concentrations in SEs, directly affecting the measured values of ionic conductivities and diffusivities. Indeed, determining the types, concentrations, dimensionalities, and effects of defects commonly present in SEs remains a great challenge in this field [64]. Computational approaches have the advantage of selectively incorporating specific defects (e.g. point defects) in SEs at specific concentrations, and monitoring their effects on the overall ion transport. We inspected the role of point defects (vacancies, interstitials, and anti-sites) on the predicted activation energies for Li + transport in Li 6 PS 5 Cl, in the context of existing experimental values.
Li + transport in Li 6 PS 5 X is mediated by intrinsic interstitials [64,71]. We showed that the existence of other defects, such as Cl S antisite and V Li provide activation energies that are also in good agreement with experimental values, suggesting that these defects may be introduced during the synthesis of Li 6 PS 5 Cl.
Plain GGA XC functionals, such as PBE and PBEsol are computationally accessible and extensively used to perform extended screening studies of SEs. Among all DFT XC functionals used during the training of MTPs, PBEsol gives results closer to experimental data for both argyrodites. We observe that the addition of optB88 vdW correction does not improve the activation energies. This result is in striking contrast with Qi et al [32], claiming that MTP trained by AIMD energy and force calculations using optB88 vdW functional are more accurate in predicting the ion transport properties for lithium superionic conductors as compared to PBE. These SEs (argyrodites) and the ones studied by Qi et al [32] (Li 3x La 2/3−x TiO 3 , Li 3 YCl 6 , Li 7 P 3 S 11 ) forms into compact structures, which typically do not require vdW functionals.
We demonstrated that changing the temperature and timescale for the AIMD simulations has direct effects on the prediction of E a values. This suggests that the simulation time and the temperature needs to be carefully chosen to get meaningful and predictive results from MTPs.
Investigating SEs with vibrational spectroscopies has, at least, two main advantages: (a) specific low-energy vibrational modes of migrating ions relate directly with their motion, and (b) they give insight into the structure of amorphous (glassy) phases that may be present in SEs. Here, we showcased that extended trajectories produced by MTPs, in particular, the autocorrelation functions can be analyzed to produce an insightful generalized density of states. This information is useful to interpret Raman, Infrared and inelastic neutron scattering experiments of SEs. GDOSs produced from extended MDs incorporate all the anharmonic effects, which are harder to capture from DFT calculations. Clearly, in the case of Li 6 PS 5 X (X = Cl, I), ensembles of modes including Li and the halide species appear at low energies.

Conclusions
In summary, we identified the main factors controlling the quality of MLPs based on the MTP when applied to the ion transport properties in SEs. Our analysis is centered on three topical SEs, Li 6 PS 5 Cl, Li 6 PS 5 I, and α-Na 3 PS 4 .
These results underline the importance of high-quality training sets in fitting MTPs and call attention to shorttimescale and high-temperature AIMD simulations to predict the room-temperature properties of materials.
We have used the derived MLPs to produce the generalized phonon density of states and van Hove correlation maps of Li 6 PS 5 Cl and Li 6 PS 5 I, which allows us to establish a direct link between the host structure and ion dynamics in these materials.