Accurate Description of Ion Migration in Solid-State Ion Conductors from Machine-Learning Molecular Dynamics

Solid-state ion conductors (SSICs) have emerged as a promising material class for electrochemical storage devices and novel compounds of this kind are continuously being discovered. High-throughout approaches that enable a rapid screening among the plethora of candidate SSIC compounds have been essential in this quest. While first-principles methods are routinely exploited in this context to provide atomic-level details on ion migration mechanisms, dynamic calculations of this type are computationally expensive and limit us in the time- and length-scales accessible during the simulations. Here, we explore the potential of recently developed machine-learning force fields for predicting different ion migration mechanisms in SSICs. Specifically, we systematically investigate three classes of SSICs that all exhibit complex ion dynamics including vibrational anharmonicities: AgI, a strongly disordered Ag$^+$ conductor; Na$_3$SbS$_4$, a Na$^+$ vacancy conductor; and Li$_{10}$GeP$_2$S$_{12}$, which features concerted Li$^+$ migration. Through systematic comparison with \textit{ab initio} molecular dynamics data, we demonstrate that machine-learning molecular dynamics provides very accurate predictions of the structural and vibrational properties including the complex anharmonic dynamics in these SSICs. The \textit{ab initio} accuracy of machine-learning molecular dynamics simulations at relatively low computational cost open a promising path toward the rapid design of novel SSICs.


I. INTRODUCTION
][10] Recent advances have increased the ion conductivities to values comparable to the ones of liquid electrolytes, making them a feasible alternative for future battery technologies. 11,12The synergy among theoretical and experimental materials studies has led to a continuous improvement of our understanding of ion-conduction phenomena, thereby generating a plethora of novel material candidates for usage of SSICs in battery devices.
4][15][16][17] Probing the ion dynamics within the solid host lattice experimentally can be achieved via impedance spectroscopy 18,19 and scattering techniques, e.g.Raman 20,21 or inelastic neutron scattering. 22,23The former technique provides important insights on the frequency range at which different ionic species migrate as well as on their densities.The latter ones complement this in important ways by providing the overall vibrational spectrum of the material, including the features associated with mobile cations and the host lattice.The interpretation of such spectra, however, is not always straightforward and often relies on suitable reference materials.Consequently, experimental efforts alone can fall short in providing a complete microscopic picture of the ion conduction mechanism.
Therefore, atomistic simulations are essential in complementing experimental investigations as they offer valuable insights into microscopic details of the conduction mechanism of mobile cations through the complex host lattice.Furthermore, they allow the virtual investigation of novel material candidates and structural phases to guide the experimental design of SSICs in a synergistic approach.Ab initio molecular dynamics (AIMD) simulations have been widely used to study the ion dynamics and particularly the coupling between mobile ions and the host lattice in SSICs, [24][25][26][27][28][29] unraveling a variety of interesting phenomena that drive the conduction of mobile cations.[40][41] Despite being a powerful and useful tool to elucidate the dynamical properties of SSICs, it is well-known that AIMD simulations suffer from a tremendous computational cost associated with the need to perform density functional theory (DFT) calculations at each time step.This limits accessible simulation times to the picosecond range and further restricts us to a small number of material systems and properties at small length scales.The prospect of machine-learning molecular dynamics (MLMD) is to overcome these shortcomings via much faster but equally accurate calculations of dynamical properties of materials. 42,43More rapid computational schemes will not only accelerate the in silico design of novel SSICs but in fact enable us to conduct more realistic calculations, e.g., regarding doped supercells.We anticipate that this will bring computational studies closer to laboratory realities and enhance their synergy with experimental synthesis and characterization. 44,45roughout an AIMD simulation, an enormous set of ab initio forces and energies is collected for atoms that typically move within a relatively narrow phase space around their equilibrium positions.This scenario appears to be ideal for machine-learning force fields (MLFFs), i.e., predicting the potential energy surface around individual atoms via supervised machine-learning (ML) techniques.9][60][61][62][63] The common target of several studies is to investigate ion-diffusion phenomena on large timescales and to exceed the structural complexity beyond classical bulk systems.Notable examples include MLFFs derived from a self-learning and adaptive database approach, which predicted anisotropic motion of Li ions in Li 10 GeP 2 S 12 (LGPS). 58Xu et al. employed moment tensor potentials for Li conductors and identified spontaneous polyanion rotation over the timescale of µs, which were however found to hamper Li-ion migration. 63Ref. 61 showed activation of PS 4 flipping events to drive structural transitions to the highly conductive Li 3 PS 4 phase at high temperatures, ruling out paddle-wheel effects in Li-ion conduction via a two-level ML scheme and ns-long MLMD simulations.Staacke et al. developed a Gaussian approximation potential to tackle ion diffusion within a lithium thiophosphate glass. 60These studies exemplify the extent by which structural complexity and/or simulation times of AIMD can be exceeded when using MLMD methods, which is important for obtaining mechanistic insights for SSICs.
Despite these incredible recent efforts, the question of whether and to what extent MLMD simulations can capture the physical mechanism behind ion transport within the strongly anharmonic potential energy surfaces of SSICs remains unaddressed.In particular, when high-throughput screening techniques shall be able to discern materials showing different ion-transport mechanisms, the question of whether MLFFs are accurate across the variety of ion-transport phenomena known to occur in SSICs is very pertinent.Here, an applicability of MLMD can deliver improved microscopic understandings of structure-property relations and the ion-migration mechanism, which will impact the design of novel material candidates.However, such prospects demand an accurate MLMD characterization of both the cation and the host lattice dynamics in SSICs.
In this work, we systematically investigate the dynamical properties of three characteristic SSICs candidates representing different transport mechanisms: AgI, a strongly disordered SSIC showing liquid-like Ag-conduction; Li 10 GeP 2 S 12 , showing concerted Li + migration; and Na 3 SbS 4 , a Na + vacancy conductor.These materials represent a variety of challenges for MLFFs including anharmonicities, highly dynamical potential energy surfaces as well as point defects.MLMD simulations were performed and compared with AIMD reference results in terms of the structural and vibrational properties of the SSIC host lattice as well as the diffusion mechanism of the mobile cations.We observe an accurate prediction of the vibrational density of states and diffusion coefficients for all systems.Moreover, even anharmonic distributions in the tilting angles within the host lattice are precisely captured.Most significantly, we demonstrate that MLMD can accurately predict the conduction mechanism of the different systems even in the presence of such extreme lattice anharmonicities.Overall, our findings demonstrate the high suitability of MLMD simulations for broad investigations of SSICs.

Ab initio molecular dynamics (AIMD) and machine-learning molecular dynamics (MLMD)
simulations were carried out with the Vienna Ab initio Simulation Package (VASP). 64,65The projector augmented wave method 66,67 and the Perdew-Burke-Ernzerhof (PBE) exchangecorrelation functional were used for AIMD simulations. 68The Brillouin zone was sampled at the Γ-point only.A canonical (N V T ) ensemble was applied in each simulation.In this work, we train MLFFs 'on-the-fly' using the machine-learning (ML) method implemented in VASP. 52,56,69After convergence of the force field, MLMD simulations are performed using the generated MLFFs to investigate their suitability for the prediction of the atomic dynamics in SSICs.The descriptors of the utilized ML method are based on a non-linear mixing of two-and three-body atomic functions describing the radial and angular distributions, respectively, within the local environment.The implemented method is a kernel-based ML model that uses a Bayesian linear regression to train on the DFT-calculated energy, forces, and stress tensors.The MLMD then makes use of a weighted sum of Kernel similarity measures, with weights obtained from the Bayesian linear regression model, to yield the prediction of forces.A summary of the ML method is given in the Supporting Information; for further details, we refer the reader to the original literature. 52,56,69Systemspecific parameters for all investigated SSICs (see Fig. 1) used in the AIMD and MLMD simulations are summarized in the Supporting Information.
We analyze the AIMD and MLFF data in two ways: first, we report the accuracy of the derived MLFFs using the Bayesian estimated error of the force (BEEF), as implemented in VASP, 52,56,69 which gives us insights into the expected forces and allows to check for convergence of the MLFFs.Additionally, forces from MLFFs are explicitly compared with DFT single-point calculations for randomly selected snapshots along the AIMD trajectories in order to quantify deviations with the root mean square error (RMSE).Hereby, we assess the quality of the MLFFs to describe the forces acting on each element separately.
Second, we calculate a variety of dynamic structural observables that are relevant for SSICs and ion conduction, starting with the atomic pair correlation function, g(r), which quantifies the short-and long-range correlation between two atomic species A and B defined as where N A and N B are the total number of A and B atoms in the supercell.Given the shortrange nature of the descriptors used in the MLFF method that is applied here, an accurate prediction of long-range behavior is an important criterion for their applicability to SSICs and will be assessed via g(r).
The vibrational density of states (VDOS) quantifies the number of modes per unit frequency range and is relevant for ion conduction because, e.g., the zero-frequency component is related to the diffusion constant, D, 70 and a presence of low-frequency vibrational modes was discussed to imply higher ion conductivity. 71,72It is computed as the Fourier transform of the normalized velocity autocorrelation function (VACF): with N and V being the number of contributing atoms and the simulation volume, respectively, and the VACF being defined as with initial time τ and v i (t) being the velocity of atom i at time t.All VDOS spectra are smoothened using a Savitzky-Golay filter with a window of 15 cm −1 for improved visibility.
The band center of a given atomic species i, ω avg,i , is calculated as the average VDOS of species i, indicated as a subscript, given as For AgI, previous work identified that significant angular fluctuations of AgI 4 units are involved in anharmonic contributions to the structural dynamics that impact ion conductivity. 26LGPS and W-doped Na 3 SbS 4 host lattices are also composed of tetrahedral units AS 4 (A = Ge, P, Sb, W), whose internal dynamics may impact the motion of mobile cations.Therefore, we calculate tetrahedral angles in AgI 4 (AS 4 ) units representing I-I-I (S-S-S) bond angles, θ i , via where d ij are atomic distances between iodide (sulfur) ions i and j within a tetrahedron.
The angular velocity of tetrahedron i is approximated as where dt is the MD simulation time step.The time evolutions of θ and ω were smoothed using a Fourier filter with a threshold frequency of 5 cm −1 to retain only the low-frequency dynamics.
Focusing on the mobile ionic species and diffusion properties, we obtain the mean-squared displacement (MSD) according to where ⟨• • • ⟩ i denotes the average over all cations i, and r i (t) gives the position of cation i at time t.D was then obtained from the MSD at final time t max as Finally, we calculate the van Hove correlation 73 and exploit it to quantify temporal correlations between cation dynamics.Specifically, the distinctive part, G d , describes the radial distribution of cations (indexed with j) with respect to a reference cation of the same (indexed with i) as a function of a time interval ∆t: with the Dirac delta function δ, the number of cations N , and the time average Conversely, the self part of the van Hove correlation, G s , provides temporal correlations of the radial distribution of cation i with itself as a function of a time interval ∆t: The relatively low computational costs of MLMD calculations allowed us to conduct five independent runs for each system in order to quantify statistical deviations of selected observables.

A. Silver iodide (AgI)
We start our analysis with silver iodide in its α-phase, AgI, which is a prototype solidstate superionic conductor: at temperatures above 420 K, AgI forms a superionic cubic phase with excellent Ag + conduction.I − ions form a body-centered cubic (bcc) lattice with a large number of interstices (6 octahedral, 12 tetrahedral, and 24 trigonal), allowing a facile motion of Ag + ions. 31,74The multitude of available Ag + sites and their strong diffusivity induce a static and dynamic disorder, and the harmonic phonon approximation is breaking down in AgI with important ramification for its ion conductivity. 26,27,75All of these render AgI an ideal, challenging system to understand the accuracy and limitations of MLFF to predict the properties and elucidate the conduction mechanisms of SSICs.
We first investigate the MLFF accuracy via the Bayesian error estimate of forces function, g(r), which shows a distinct I-Ag peak at 2.8 Å and a significant broadening at larger distances, see Fig. 2b.I − ions exhibit long-range order, while the Ag-Ag g(r) shows a liquid-like behavior with a value of ∼1 across the whole range.Remarkably, MLMD simulations reproduce the AIMD-calculated g(r) of both the iodide host lattice and the disordered Ag + ions.In view of the very different dynamic characteristics of the mobile Ag cation and I host lattice, this finding is important and encouraging.We further investigate the VDOS (see Fig. 2c and Methods section for computational details) and find two broad peaks around 25 and 100 cm −1 .Contributions from Ag + and I − ions dominate the former, while the latter is dominated by I − modes, in good agreement with previous reports. 26,76We further obtain non-zero values at 0 cm −1 due to slow Ag + vibrational modes, which is in line with the high ion conductivity of AgI.Our MLMD simulations show an excellent agreement over the whole frequency range from 0 to 200 cm −1 , suggesting that MLFFs can accurately capture and reproduce not only structural but also vibrational properties even in the highly complex potential energy surface of AgI.
We now focus on the Ag + ion conduction mechanism and investigate the MSD of Ag + ions, see Fig. 2d and Methods section for computational details.We observe a close to linear relationship for MSD(t) in both AIMD and MLMD simulations.Notably, the linear increase of MSD with time is known from Brownian motion, 77 which together with g(r) ∼ 1 suggests random dynamics of Ag + ions in a liquid-like environment.D Ag + values extracted from the MSD from AIMD and MLMD simulations are 2.61 × 10 −5 cm 2 s −1 and 2.99 × 10 −5 cm 2 s −1 , respectively, which are in reasonable agreement and in line with previous studies. 78,79cause anharmonic relaxational dynamics of the iodide host-lattice were discussed as a driving force for Ag + conduction, 26 it is important to examine whether MLFFs can capture such dynamic effects.First, in both AIMD and MLMD we observe large deviations of I-I-I bond angles, θ, from the expected value of 70.5 • to below 45 • and above 100 • , see Fig. 3a.As shown by Brenner et al., 26 these distortions remain present for long timescales, as also seen in the slow relaxational components of the θ angles we compute (see Fig. 3a).Interestingly, AIMD and MLMD simulations show a comparable extent of deviations from the expected value.A statistical analysis of the angle probability density, ρ(θ), taken from 48 different tetrahedral angles in the supercell along the trajectory, shows that values between 30 • and 110 • occur in the simulations (see Fig. 3b).ρ(θ) as computed in AIMD shows a more skewed form, while the counterpart from MLMD is slightly more pronounced at 70.5 • .The tail at small values of θ is captured well by MLMD but a more pronounced disagreement is observed for larger values of θ.The standard deviation among several MLMD runs is found to be small, which indicates a high degree of reproducibility of these simulations for capturing such distinct local distortions.
Because of the highly dynamic nature of AgI, the restoring force acting on iodide ions are anticipated to be rather small, i.e. small ω values may be predicted.This is in line with the results shown in Fig. 3a, where deviations in θ from the expected value occur for long timescales of ∼ 20 ps.Indeed, we find small values of ω between −0.2 • fs −1 and 0.2 • fs −1 (see Fig. 3c), confirming the weak restoring forces.MLMD simulations accurately capture the AIMD probability density of angular velocities, ρ(ω), over the entire range suggesting an accurate prediction of the dynamical changes in the iodide host lattice (see Fig. 3c).
A potential reason for the deviations in the θ distribution of MLMD compared to AIMD may be rooted in the limited phase space that is being sampled during MLFF training.To investigate this hypothesis, we retrained the MLFF at an increased temperature of 800 K and performed another MLMD simulation at 500 K using the alternative MLFF.In this way, the system explores a larger part of phase space during MLFF training, expanding into strongly anharmonic regions due to the increased thermal fluctuations.However, the deviation in ρ(θ) is found to be similar even when training is conducted at higher temperatures (Figure S3, Supporting Information).Consequently, we may assign the remaining deviations to the choice of two-and three-body descriptors in the MLFFs, for which more complex ones may be required to capture all aspects of ρ(θ) with a high accuracy.
Finally, we investigate the van Hove correlation function of Ag + ions, considering the self part, G s and the distinctive part, G d , see Fig. 4 and Methods section for the computational details.In line with the monotonous increase of the MSD (cf.Fig. 2d), an increasing width Lithium-based SSICs emerged as a promising alternative to replace their liquid counterparts without changing the full battery stack because of their excellent ion conductivities. 802][83][84][85][86][87][88][89][90] Many studies have investigated the underlying mechanism that drives Li + migration through the solid host lattice, 87,[91][92][93] and AIMD studies suggested a concerted migration of Li + ions to be at play. 39,94Specifically, the coupled motion of multiple Li + ions can lower migration energy barriers significantly compared to a scenario where ions migrate individually.This mechanism of coupled Li + motion is similar to the situation in liquid electrolytes 95,96  in particular for Li-Li and Li-S.However, we do observe differences in the Li-P and Li-Ge g(r) peak intensities, whereas their peak positions are accurately captured.Specifically, their first and second g(r) peaks are not precisely reproduced by MLMD.These deviations are in line with the enhanced RMSE values for Ge and P atoms, which is likely caused by the small amount of Ge and P atoms in the supercell, i.e. 8 and 16 atoms, respectively, and consequently the low number of local reference configurations for each type, as mentioned in the computational details.It is noted that MLMD further captures the g(r) within the polyanions of the host lattice (see Figure S5, Supporting Information).
We now investigate the VDOS of LGPS (see Fig. 5c) and obtain good agreement between AIMD and MLMD simulations for the entire frequency range.All element-wise contribu-tions, particularly the ones of Li + ions, are closely mirrored in MLMD.Li vibrations appear throughout the whole frequency range up to 650 cm −1 , dominating in particular the region between 200 cm −1 to 400 cm −1 , whereas contributions by Ge and P are found mostly for higher-frequency modes at 550 cm −1 .There are some deviations between AIMD and MLMD visible in the high-frequency range of the polyanion data, which are likely not very relevant because ion migration is driven by lower-frequency modes.Furthermore, we observe a finite VDOS value at zero frequency, pointing to the Li + migration.MLMD simulations, averaged across 5 runs, show a continuous increase in the MSD (Fig. 5d).In contrast to our findings for AgI, we notice that the variation in the Li + MSD is substantial, which suggests an increased stochastic nature of the ion-migration mechanism.The MSD from the single AIMD run follows the MLMD one to a good extent.D Li + values of 5.33 × 10 −5 cm 2 s −1 and 4.27 × 10 −5 cm 2 s −1 from AIMD and MLMD simulations, respectively, were extracted, which are in line with previous computational studies. 97We assign the observed deviations to the increased stochastic character of Li migration in LGPS.Thus, large simulation times -exceeding feasible AIMD timescales -are required to predict reliable MSD values.
We now analyze the migration mechanism of Li + ions.Here, we calculate the self and distinctive part of the van Hove correlation function for Li + ions in LGPS from AIMD and MLMD trajectories.MLMD results of G d (r, ∆t) show a substantial correlation at 0 Å starting at ∼ 1 ps.This means that a Li + ion follows a migrating one within timescales of around 1 ps, supporting the concerted motion of Li + ions (see Fig. 6 and Figure S6, Supporting Information).Strong correlations exist around 4 Å and 6.5 Å, confirming that the long-range order of Li + ions remains unaffected even though Li + ions migrate, in line with previous results. 39G s (r, ∆t) shows a continuous increase in correlation during the simulations, showing that the Li + ions are migrating, which is consistent with the MSD (see Fig. 5d).In contrast to AgI, we observe a substantial correlation at r = 0 Å for all ∆t values.This underpins that a substantial amount of Li + ions remain located near their equilibrium positions, while all Ag + ions contribute to the overall conduction in AgI.Notably, the MLFF shows a fair agreement for G s (r, ∆t) and G d (r, ∆t) with the AIMD reference (see Figure S6, Supporting Information).Deviations in the overall values and intensities are caused by the limited available timescales from AIMD simulations.Importantly, we can still confirm that MLMD captures the concerted migration mechanism in LGPS.
We now focus on the dynamics of the PS 4 and GeS 4 tetrahedra in the LGPS host lattice because distortions in the polyanion backbone are likely to play a crucial role in facilitating cation diffusion. 98Moreover, it is crucial to investigate the aforementioned consequences of the increased RMSE of the MLFF for Ge and P, as well as the deviations seen in g(r), for the dynamical properties of the host lattice.We compute the distribution of θ angles within the PS 4 and GeS 4 tetrahedra and compare the AIMD and MLMD probability densities, ρ(θ), see Fig. 7.The agreement between AIMD and MLMD is overall good.We generally find small variations in the θ angles suggesting relatively rigid tetrahedra.PS 4 tetrahedra show fluctuations that are largely harmonic in the range between 50 • and 70 • .Interestingly, the angles in GeS 4 tetrahedra exhibit a distribution that is significantly more asymmetric with an extended tail toward large values.This suggests the presence of substantial anharmonicities in the rotational dynamics of GeS 4 tetrahedra.The angular velocity distribution, ρ(ω), further underlines the accuracy of the MLMD simulations in capturing not only the angular distribution but also their dynamical changes (Fig. 7).The higher RMSE values of the forces from MLFFs for P and Ge atoms (see Figure S4, Supporting Information) may explain the remaining deviations present in the angular distributions.In contrast, the RMSE of Li forces is considerably lower, in line with our finding that Li + migration is captured well.
Having assessed the dynamical properties of mobile cations and host lattice separately, we now investigate the interaction between PS 4 and GeS 4 polyanions and adjacent Li + ions.
An accurate description from MLMD is relevant for capturing pertinent dynamic coupling phenomena, such as paddle-wheel effects, which potentially drive Li + ion migration. 30,34,63 compute the 2D probability distribution, ρ(r Li-S , ϕ), which measures occurrences of Li-S distances, r Li-S , as a function of P-S-Li bond angles, ϕ, see also Figure S7 Na-based SSICs gained significant attention as a potential replacement for their Libased counterparts.Sodium thiophosphates, Na 3 PS 4 , 99,100 and sodium thioantimoniates, Na 3 SbS 4 , [101][102][103] have been the focus of many studies in this context because of their integrability into all-solid-state sodium batteries and due to their high ion conductivities. 104,105reover, a variety of interesting findings have been reported on the atomic dynamics in these materials.This includes vibrational anharmonicity, 28,36 a rich polymorphism, 24,106 and coupled dynamics of anions and cations that have been discussed in the context of paddlewheel effects. 24,32A peculiarity of current Na-based SSICs is their tunable ion conductivity by introducing Na vacancies via aliovalent doping.Several studies suggested that Na + diffusion in pristine materials remains absent even at temperatures as high as 900 K and only becomes viable upon introduction of defects. 107,108Such can be achieved through a variety of different dopants, including halides [109][110][111] and tungsten, 103,112,113 which commonly substitute a pnictogen or chalcogen, respectively, and leave Na + vacancies behind to counterbalance the change in oxidation state.Na + ions may then diffuse via vacancy migration with DFTcomputed migration barriers as low as 0.1 eV, resulting in ion conductivities in the mS cm −1 range. 114MD has been instrumental in demonstrating the migration of Na + vacancies within the network of [PnCh 4 ] 3− (Pn={P, Sb}; Ch={S, Se}) polyanions. 24,28,108,111Specifically, the anharmonic character of the host-lattice vibrations and local symmetry breaking, which cannot be captured by harmonic phonon calculations alone, have motivated use of finitetemperature MD simulations. 36,71,106,115However, computational costs associated with doped scenarios of, e.g., Na 3 PnCh 4 SSICs, are very high.First, large supercells are needed to mimic small doping concentrations at the Pn site.Second, Na + migration occurs by hopping across vacancies, which are rare events despite modest ion conductivities and low migration barriers found in these materials.Consequently, long simulation times are needed to sample sufficiently many hopping events and describe the coupling between host lattice and migrating Na + ions.
MLMD trained via AIMD data may well constitute a promising way to overcome these computational challenges, potentially allowing for the study of extended length and time scales, which we investigate here for the scenario of aliovalent doping.Specifically, we focus on tungsten-doped Na 3−x Sb 1−x W x S 4 , which was shown to attain excellent ion conductivities. 112,113The mechanism involves W substituting a Sb site, forming [WS 4 ] 2− polyanions that exhibit different charges than the original [SbS 4 ] 3− ones, and a variation in the chemical environment of S atoms.As a consequence, Na vacancies, V − Na , are introduced which may diffuse within the crystal.The local variation in chemical bonds and charge as well as changes in the potential energy surface due to V − Na migration, can be expected to result in challenges for the training procedure of MLFFs.
We start by considering the accuracy of MLMD simulations for pristine Na 3 SbS 4 at T = 300 K.The MLFF quickly converges within 1000 steps, showing maximum force errors of below 0.02 eV Å−1 after 200 steps already (see Figure S9, Supporting Information).This results in precise predictions of g(r) among all elements in the anion host lattice and the Na + ions (see Figure S10, Supporting Information).The VDOS shows the expected polyanion vibrations at larger frequencies of 350 cm −1 to 400 cm −1 , and broad contributions in the range of ≤ 200 cm −1 caused in particular by Na and S vibration (see Figure S11, Supporting Information).The latter suggests the presence of anion-cation coupling, in line with previous studies. 28,36,106Notably, the MLMD data capture mode frequencies across the entire frequency range of vibrations, but moderate differences in their intensities remain present.
We further observe that Na + ions merely oscillate around their lowest-energy positions, while migration throughout the pristine lattice remains absent in our AIMD and MLMD simulations (see Figures S12-S14, Supporting Information), as previously reported. 107xt, we investigate the accuracy of MLFFs for the tungsten-doped Na 2.94 Sb 0.94 W 0.06 S 4 .
A new set of force fields is trained starting from the MLFF of the pristine Na 3 SbS 4 .The BEEF values during the training procedure largely remain below 0.02 eV Å−1 after around 1000 steps (see Figure S15, Supporting Information).As expected, the error of the MLFF is higher than the one for the pristine system due to increased variations in the local co- An important question is whether MLMD simulations accurately describe the vibrational properties of W-doped Na 2.94 Sb 0.94 W 0.06 S 4 .Indeed, we observe good agreement between MLMD and AIMD data for the overall shape of the VDOS and individual elementwise contributions in the region ≤ 300 cm −1 (see Fig. 8b) despite substantial changes in the local chemical environments.However, the intensity of higher-energy sulfide modes shows larger deviations despite their frequencies being captured well, as explicitly shown in ∆VDOS plots (Figure S18, Supporting Information).While an accurate representation of low-frequency modes is more important for describing SSICs, reducing such deviations for the high-frequency vibrational regime could be achieved with improved MLFFs.We also investigate the Na band center, given by the average vibrational frequencies of the projected Na VDOS (see Eq. ( 4)), because it was recently identified as a good indicator for changes in the overall lattice stiffness. 71Here, AIMD and MLMD simulations give Na band center values of ω avg,Na = 127.5 cm −1 and 126.5 cm −1 , respectively, showing excellent agreement.
By contrast, the pristine Na 3 SbS 4 exhibits a larger value of 129.7 cm −1 (129.1 cm −1 ) from MLMD (AIMD) simulations, which shows that a small softening of the lattice occurs upon introducing W dopants and that this correlates with the higher Na + conductivity.
The prediction of Na + diffusion coefficients from AIMD simulations is greatly complicated by random Na + hopping processes.Consequently, long simulations and/or large cell sizes are required to obtain reliable MSD values, which again underlines the need for reliable MLMD.Our analysis of the Na + MSD from MLMD simulations highlights these complications by showing a wide spread of MSD values over time (Fig. 9a), from which we obtain D + Na = (1.13 ± 0.69) × 10 −6 cm 2 s −1 .The MSD computed from the single AIMD trajectory lies within the standard deviation of MLMD values.Van Hove correlations for Na 2.94 Sb 0.94 W 0.06 S 4 exhibit a behavior that significantly deviates from AgI and LGPS: G s develops an additional region of strong correlations after 7 ps at approximately 3.5 Å, corresponding to the distance between neighboring Na + sites and confirming the migration of Na + ions (see Fig. 9b and Figure S19, Supporting Information).Notably, the region between the two strongly correlated G s features is only weakly populated, which points to the migration of Na + ions via thermally activated hopping.Concerted migration can be ruled out from G d (see Fig. 9c), showing only a slowly increasing component after 6 ps.
Because several studies suggested that coupling between cations and host lattice is a key factor in facilitating cation diffusion, 24,32 we investigate the dynamics of SbS  MLMD data were averaged across 5 independent trajectories.anharmonicity in the rotational motions of SbS 4 tetrahedra.This is also visible in the angular distribution, ρ(ω), of SbS 4 tetrahedra, which shows small deviations around −0.5 and 0.5 • fs −1 from harmonic behavior.Notably, ρ(ω) of WS 4 tetrahedra shows a more harmonic distribution.The overall small deviations between AIMD and MLMD are found to be more pronounced for WS 4 , which likely is a result of the small sampling size of WS 4 .
We now focus on correlations between SbS 4 and WS 4 tetrahedra and neighboring Na + ions.We obtain the 2D probability density ρ(r Na−S , ϕ), by defining the distance, r Na−S , between neighboring Na + ions and S atoms, and the dihedral angle ϕ, see Fig. 10   given more space to migrate in the proximity of the W-dopant.The absence of correlated branches as we found in case of LGPS (see Fig. 7) can be rationalized by the nature of Na + migration discussed above when investigating G s for Na 2.94 Sb 0.94 W 0.06 S 4 (see Fig. 9): Na + cations spend comparably little time in-between equilibrium sites, which is characteristic of a thermally-activated cation hopping, whereas Li + ions in LGPS spend more time in-between their equilibrium sites.Altogether, these findings show that the host lattice dynamics of Na 2.94 Sb 0.94 W 0.06 S 4 are correlated with cation diffusion and suggest that ion conduction can be tuned through adjusting the dynamics of the host lattice.

IV. DISCUSSION
Our systematic analysis found that MLMD simulations are accurate for the vibrational properties of mobile cations and the host lattice.Anharmonic properties are predicted well, and the interactions between mobile cations and the polyanionic backbone are precisely captured for all investigated systems.These results are especially intriguing as various ionconduction mechanisms were investigated: liquid-like behavior of Ag + ions in AgI, concerted migration of Li + in LGPS, and thermally activated V − Na vacancy hopping in Na 3 SbS 4 .To the best of our knowledge, a systematic investigation of MLFF accuracy covering such diverse ion migration phenomena in SSICs has not been reported yet.The majority of existing studies rely on the accuracy of MLFFs in comparison to DFT reference calculations.However, recent reports highlight that benchmarking the accuracy of forces alone does not necessarily align with an accurate prediction of relevant material properties. 116Here, both structural and vibrational properties of mobile cations and host lattice are found to be captured well, providing confidence that the utilized MLFFs are highly suitable for a broad range of SSICs.
We emphasize that such accurate predictions from the exploited MLMD model were not necessarily expected.The locality of descriptors, which only include atoms that are located within a certain cutoff radius, is often discussed as one downside for larger systems next to advantages such as simplicity and transferability. 42,1179][120][121] In the case of SSICs, the short-range nature of interactions and the large dynamical distortions exhibited by the ions likely aids the kernel-based ML procedure.A dominance of short-range interactions can lead to long-range phenomena being less relevant effectively, which can rationalize our accurate MLMD simulation results.Correlated ion phenomena at short ranges still remain well represented as we have demonstrated in the case study of concerted migration of Li ions in LGPS.
Finally, we would like to provide a compact perspective on the relevance of our findings for future development of materials and the role MLMD simulations may play in this regard.
Fast and reliable MLMD simulations with ab initio accuracy can open novel possibilities when guiding experimental materials design and characterization.Raman spectroscopy, for example, represents one of the most widely used, non-destructive techniques for providing vibrational fingerprints of a solid-state material.Even though the first approaches towards high-throughput Raman calculations and characterization have been made recently, [122][123][124][125][126] first-principles calculations of Raman spectra remain essential for the interpretation of experimental measurements and prediction of vibrational properties in novel material systems.
Commonly used Raman calculations based on the harmonic approximation are, however, not sufficient to capture the entire vibrational dynamics of materials that are as anharmonic as the investigated SSICs here.These inherent limitations can be overcome through molecular dynamics simulations that account for the whole anharmonicity of such material, [127][128][129] which in the past found only limited attention due to the tremendous computational costs associated with AIMD.The excellent accuracy of the presented MLMD results for the variety of vibrational properties in the investigated SSICs may render calculations of finitetemperature Raman spectra viable and hereby inspire novel data-driven approaches, combining ML-models for Raman activities [130][131][132][133][134][135][136] with MLMD to launch a flourishing synergy between experiments and computation.
1][142][143] Incorporating vibrational properties as well as information about the ionconduction mechanism derived from MLMD simulations may improve the selectivity of existing screening approaches and potentially deepen our understanding of the structureproperty relations of SSICs.

V. CONCLUSION
In summary, we have provided a systematic investigation of the accuracy of MLMD simulations to predict the dynamical properties associated with ion conduction across various mechanisms in a broad range of SSICs.Results from MLMD were compared to AIMD reference data using several measures that include the actual force errors, structural and dynamical properties of mobile cations as well as of the host lattice, and correlations among the motion of cations and rotational degrees of freedom in the polyanion backbone.Our results demonstrate that precise predictions of dynamical properties of SSICs exhibiting different transport mechanisms that involve liquid-like cation motion in AgI, concerted cation migration in LGPS, and thermally-activated vacancy hopping in W-doped Na 3 SbS 4 are feasible.Through our MLMD simulations, we could further provide a detailed analysis of the anharmonic character of host lattice vibrations, which take a significant role in enabling a facile migration of mobile cations.Our findings are very encouraging and may further enhance trust in MLMD simulation results and their suitability for a broad range of SSICs with different physical mechanisms governing ion migration and structural variability.We believe that the increasing availability of MLMD methods can enhance synergies among experiments, computations, and material design efforts.An exciting avenue along this quest revolves around the development of novel data-driven screening techniques for SSICs, which will be aided by finite-temperature materials predictions through MLMD simulations at relatively low computational cost.

FIG. 1 .
FIG. 1. Schematic representation of the solid-state ionic conductors investigated in this work: (Left) AgI, with Ag in silver and I in purple; (mid) Li 10 GeP 2 S 12 (LGPS) with Li in bright orange, Ge in grey, P in green and S in yellow; (right) Na 3−x Sb 1−x W x S 4 with Na in dark orange, W in blue, Sb in purple and S in yellow.The ion-conduction mechanisms of each system are sketched in the lower panels.

(FIG. 2 .
FIG. 2. (a) Bayesian error estimate of forces (BEEF) as function of the number of training steps, N train , during training of the MLFF for AgI.Panels (b) and (c) visualize the atomic pair correlation function, g(r), and the vibrational density of states (VDOS), respectively.In panel (d) the mean square displacement (MSD) of Ag + ions is shown as a function of time, t, with the dashed line being the average of 5 MLMD runs and the blue-shaded region the standard deviation.

FIG. 3 .
FIG. 3. (a) Representative time evolution of I-I-I bond angle, θ, for selected AgI 4 tetrahedra (see Figure S2, Supporting Information, for schematic representation).Thick (thin) lines in panel (a) represent the Fourier filtered (unfiltered) values of θ using a cutoff frequency of 5 cm −1 .(b) Probability density, ρ(θ), as obtained from 48 tetrahedra within the AgI supercell.(c) Probability densities of the angular velocity, ρ(ω), calculated from the time evolution of θ.For ρ(θ) and ρ(ω), data averaged from 5 independent MLMD runs are shown and the light-blue area indicates the standard deviation.

FIG. 4 .
FIG. 4. Spatiotemporal correlations of Ag + ion migration in AgI: (a) and (c) visualize the self-part, G s (r, ∆t), of the van Hove correlation extracted from AIMD and MLMD simulations, respectively.Panels (b) and (d) visualize the distinctive part, G d (r, ∆t), from AIMD and MLMD simulations, respectively.MLMD data were averaged across 5 independent trajectories.

FIG. 5 .
FIG. 5. (a) BEEF as function of N train during training of the MLFF for LGPS.Panels (b) and (c) visualize g(r) and the VDOS, respectively.In panel (d) the MSD of Li + ions is shown as a function of t, with the dashed line being the average of 5 MLMD runs and the blue-shaded region the standard deviation.
, Supporting Information, for a graphical representation.Highest correlation between PS 4 tetrahedra and Li ions is found in a small range of 80 • ≤ ϕ ≤ 95 • at small Li-S distances of ≤ 3 Å (see Fig. 7 and Figure S8, Supporting Information).Two less intense branches are found leaving this highly correlated region.The first one reaches higher ϕ angles of up to 130 • with Li-S distances spread over the whole considered range.The second, less-intense branch simultaneously goes to smaller ϕ angles and larger Li-S distances.The correlation of Li-S distances with GeS 4 tetrahedra shows a squeezed form with small changes in the distances, showing ϕ angles of 75 • to 110 • .This suggests a larger rotational degree of freedom of Li ions around GeS 4 than around PS 4 tetrahedra, which we may assign to the larger mass and less negative charge of the [GeS 4 ] 2− polyanions.Interestingly, Li-S distances increase only when ϕ ≤ 70 • or ϕ ≥ 120 • , representing a potential indicator for the design of host lattice structures to enable Li ion migration.Direct comparison with AIMD results (see Figure S8, Supporting Information) shows an excellent match in the probability distributions with MLMD data, supporting the accurate prediction of interactions between mobile cations and the host lattice.These results underline that MLMD captures and reproduces the structural, vibrational, and concerted nature of Li ion migration in LGPS.C. Sodium Thioantimonate (Na 3 SbS 4 )
4 and WS 4 tetrahedra in the host lattice of Na 2.94 Sb 0.94 W 0.06 S 4 .Distributions of θ angles within SbS 4 and WS 4 tetrahedra are obtained in order to compare the AIMD and MLMD probability densities, ρ(θ), see Fig.10a and d, respectively.The overall agreement between AIMD and MLMD is again satisfactory.The overall small dynamic changes of θ indicate that the tetrahedra form relatively rigid structures.Notably, θ of SbS 4 tetrahedra exhibit an asymmetric distribution with an extended tail toward larger angle values whereas θ of WS 4 tetrahedra exhibit an essentially symmetric characteristic.This indicates the existence of

FIG. 9 .
FIG. 9. (a) MSD of Na + ions in Na 2.94 Sb 0.94 W 0.06 S 4 as a function of t, with the dashed line being the average of 5 MLMD runs and the blue-shaded region the standard deviation.Spatiotemporal correlations of Na + ion migration: Panels (b) and (c) visualize G s (r, ∆t) and G d (r, ∆t), respectively.
for MLMD and Figure S20, Supporting Information, for AIMD results.The largest correlation between SbS 4 tetrahedra and Na ions exists in a small range of 80 • ≤ ϕ ≤ 100 • with r Na-S remaining below ≤ 4 Å.Minor intensities are found in two regions: one in a range of 100 • ≤ ϕ ≤ 120 • and r Na−S of ≤ 4 Å, and the other in 60 • ≤ ϕ ≤ 80 • and 5 Å ≤ r Na−S ≤ 6 Å.The noticeable intensities at lower r Na−S values showing a broadened ϕ distribution between 80 • and 120 • indicate that Na + ions experience a significant degree of rotational freedom around SbS 4 tetrahedra.In regard to WS 4 tetrahedra, the intensity in the range 60 • ≤ ϕ ≤ 80 • and 5 Å ≤ r Na−S ≤ 6 Å is found to be more pronounced than in case of SbS 4ones.This observation suggests a facilitated ion migration when Na + cations are passing by WS 4 tetrahedra, which can be explained by W-S bond lengths being shorter than the Sb-S ones and the fact that WS 4 tetrahedra occupy less space.Consequently, Na + ions are