Takeru
Miyagawa
,
Namita
Krishnan
,
Manuel
Grumet
,
Christian Reverón
Baecker
,
Waldemar
Kaiser
* and
David A.
Egger
*
Department of Physics, TUM School of Natural Sciences, Technical University of Munich, 85748 Garching, Germany. E-mail: waldemar.kaiser@tum.de; david.egger@tum.de
First published on 9th April 2024
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; Na3SbS4, a Na+ vacancy conductor; and Li10GeP2S12, which features concerted Li+ migration. Through systematic comparison with 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 ab initio accuracy of machine-learning molecular dynamics simulations at relatively low computational cost opens a promising path toward the rapid design of novel SSICs.
Deriving structure–property relations between SSIC compositions and their ionic conductivities is key when establishing design rules for screening potential SSIC candidates in experimental and theoretical high-throughput approaches.13–17 Probing the ion dynamics within the solid host lattice experimentally can be achieved via impedance spectroscopy18,19 and scattering techniques, e.g. Raman20,21 or inelastic neutron scattering.22,23 The 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–29 unraveling a variety of interesting phenomena that drive the conduction of mobile cations. Such computational studies have suggested various mechanisms, including the ‘paddle-wheel effect’ that is coupling cation translational and polyanion rotational motions to impact cation conduction.30–34 Furthermore, various SSICs were found to exhibit anharmonic host lattice dynamics, which involve the mobile ions being driven into strongly anharmonic regions of the potential energy surface.26,35–37 Moreover, the collective motion of cations through the soft lattice, in the form of concerted migration, was found to substantially alter the local potential barriers experienced by cations for a variety of lithium and sodium ion conductors.38–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,43 More 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,45
Throughout 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. Various MLMD methods have therefore been developed to mitigate the intrinsic limitations of AIMD.46–52 Among them, on-the-fly active learning methods can update the MLFF by training with additional DFT calculations when facing structures of low certainty that have not been sufficiently sampled during the training steps.52–57
The enormous recent progress in the area of MLMD has already had an impact in the SSIC community, e.g., several recent studies employed MLMD to investigate ion diffusion properties and the role of the host lattice.58–64 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 Li10GeP2S12 (LGPS).58 Xu 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.62 Ref. 63 showed activation of PS4 flipping events to drive structural transitions to the highly conductive Li3PS4 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.60 These 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; Li10GeP2S12, showing concerted Li+ migration; and Na3SbS4, 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.
In this work, we train MLFFs ‘on-the-fly’ using the ML method implemented in VASP.52,56,70 After 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 ESI;† for further details, we refer the reader to the original literature.52,56,70 System-specific parameters for all investigated SSICs (see Fig. 1) used in the AIMD and MLMD simulations are summarized in the ESI.†
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,70 which gives us insights into the expected force accuracies 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
(1) |
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,71 and a presence of low-frequency vibrational modes was discussed to imply higher ion conductivity.72,73 The VDOS of an atomic species i is computed as the Fourier transform of the velocity autocorrelation:
(2) |
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
(3) |
For AgI, previous work identified that significant angular fluctuations of AgI4 units are involved in anharmonic contributions to the structural dynamics that impact ion conductivity.26 LGPS and W-doped Na3SbS4 host lattices are also composed of tetrahedral units AS4 (A = Ge, P, Sb, W), whose internal dynamics may impact the motion of mobile cations. Therefore, we calculate tetrahedral angles in AgI4 (AS4) units representing I–I–I (S–S–S) bond angles, θi, via
(4) |
(5) |
Focusing on the mobile ionic species and diffusion properties, we obtain the mean-squared displacement (MSD) according to
(6) |
(7) |
Finally, we calculate the van Hove correlation74 and exploit it to quantify spatiotemporal correlations between cation dynamics. Specifically, the distinctive part, Gd, describes the radial distribution of cations (indexed with j) with respect to a reference cation of the same type (indexed with i) as a function of a time interval Δt:
(8) |
(9) |
The relatively low computational costs of MLMD calculations, reported in Table S1, ESI,† allowed us to conduct five independent runs for each system in order to quantify statistical deviations of selected observables.
We first investigate the MLFF accuracy via the Bayesian error estimate of forces (BEEF) during the training period (Fig. 2a), showing well-converged force predictions to ≤0.05 eV Å−1. The actual root mean squared errors (RMSE) of the forces in the MLFF, referred to explicit DFT force calculations, are 45.5 meV Å−1 and 40.2 meV Å−1 for I and Ag atoms, respectively, see Fig. S1, ESI.† Next, we consider the AIMD-computed (see Methods section for computational details) atomic pair correlation 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,77 We 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,78 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 cm2 s−1 and 2.99 × 10−5 cm2 s−1, respectively, which are in reasonable agreement and in line with previous experimental and theoretical studies.79,80
Because 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.
Fig. 3 (a) Representative time evolution of I–I–I bond angle, θ, for selected AgI4 tetrahedra (see Fig. S2, ESI,† 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. |
To further examine the θ distributions, we compute the angular velocity, ω, see eqn (5). 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. This is a known strategy to reduce the risk of excessive extrapolation during MLMD simulations.42,81 However, the deviation in ρ(θ) is found to be similar even when training is conducted at higher temperatures (Fig. S3, ESI†). 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.
We now investigate the van Hove correlation function of Ag+ ions, considering the self part, Gs and the distinctive part, Gd, 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 in the Gs values with time is observed. Interestingly, Gs values at r = 0 Å quickly drop to zero as Δt is increased. This implies that Ag+ ions barely return to their previous location, which is a further aspect confirming the randomness in the Ag+ ion dynamics. Furthermore, after fast broadening of the strongly correlated Gs domain within Δt ≤ 0.5 ps, its center increases continuously with time. The diffusive character of Ag+ dynamics is visible as a rapid and continuous increase of Gs values that reflect traveled distances of Ag+ ions, contributing to the remarkable ion conductivity of AgI. At the same time, no significant Gd peaks are observed. Rather, the Gd data confirm the fully homogeneous distribution of cations over space and time without significant temporal correlations between different Ag+ ions. Importantly, MLMD simulations accurately capture the shape and intensities of Gs and Gd as obtained in AIMD. Overall, our results suggest that MLMD is able to capture the structural, vibrational, diffusive, and angular properties of AgI even with the intricate potential energy surface this material experiences at higher temperatures.
Finally, we investigate size transferability and apply the generated MLFF from the 128 atom supercell to predict the structural and dynamical properties of a 1024 atom supercell of AgI. We find an excellent match in the g(r), VDOS, and Ag+ MSD, see Fig. S4.† Notably, we observe that the Ag+ MSD gets more linear, which is caused by the larger amount of Ag+ ions that can move within the larger supercell. The findings suggest a transferability of MLFFs to larger systems sizes for cases where the crystal structure remains unchanged.
As above, we analyze the accuracy of the trained MLFF first via the force errors during the training period (Fig. 2a) and observe fast convergence of the BEEF to values below 0.05 eV Å−1 (see Fig. 5a). For the individual chemical species, we observe RMSE values of MLFF with respect to DFT force values of 27.7 meV Å−1 for Li, 75.7 meV Å−1 for Ge, 85.0 meV Å−1 for P, and 45.4 meV Å−1 for S (see Fig. S5, ESI†). The calculated g(r) show a first-order peak around 2.5 Å for Li–S and broadened higher-order peaks (see Fig. 5b). The g(r) data for Li–P and Li–Ge both show a peak at 3.5 Å and exhibit an overall structure that is consistent with the expected long-range order of the host lattice. By contrast, a suppressed long-range order is found for Li–Li with g(r) ≈ 1 across the whole range. MLMD simulations accurately capture the features of g(r) for Li with most elements, in particular for Li–Li and Li–S. However, we do observe small differences in the Li–P and Li–Ge g(r) peak intensities, whereas their peak positions are accurately captured. Specifically, their second g(r) peaks are not precisely reproduced by MLMD. These deviations are in line with the slightly 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 Fig. S6, ESI†).
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 contributions, particularly the ones of Li+ ions, are closely mirrored in MLMD. Li vibrations appear throughout the whole frequency range up to 650 cm−1, being strongly coupled to S vibrations in particular in 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−6 cm2 s−1 and 4.50 × 10−6 cm2 s−1 from AIMD and MLMD simulations, respectively, were extracted, which are in line with previous computational studies99 and quasi-elastic neutron scattering measurements.100 We 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 Gd(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 S7, ESI†). 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.39Gs(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 MLMD shows a fair agreement for Gs(r, Δt) and Gd(r, Δt) with the AIMD reference (see Fig. S7, ESI†). 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.
Fig. 6 Spatiotemporal correlations of Li+ ion migration in LGPS: Panels (a) and (b) visualize Gs(r, Δt) and Gd(r, Δt), respectively. MLMD data were averaged across 5 independent trajectories. |
We now focus on the dynamics of the PS4 and GeS4 tetrahedra in the LGPS host lattice because distortions in the polyanion backbone are likely to play a crucial role in facilitating cation diffusion.101 Moreover, it is interesting 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 PS4 and GeS4 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. PS4 tetrahedra show fluctuations that are largely harmonic in the range between 50° and 70°. Interestingly, the angles in GeS4 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 GeS4 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 Fig. S5, ESI†) 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 PS4 and GeS4 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,62 We compute the 2D probability distribution, ρ(rLi–S, ϕ), which measures occurrences of Li–S distances, rLi–S, as a function of P–S–Li bond angles, ϕ, see also Fig. S8, ESI,† for a graphical representation. The highest correlation between PS4 tetrahedra and Li ions is found in a small range of 80° ≤ ϕ ≤ 95° at small Li–S distances of ≤3 Å (see Fig. 7 and S9, ESI†). 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 GeS4 tetrahedra exhibits 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 GeS4 than around PS4 tetrahedra, which we may assign to the larger mass and less negative charge of the [GeS4]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 Fig. S9, ESI†) 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 and vibrational properties, as well as the concerted nature of Li ion migration in LGPS.
AIMD has been instrumental in demonstrating the migration of Na+ vacancies within the network of [PnCh4]3− (Pn = {P, Sb}; Ch = {S, Se}) polyanions.24,28,111,114 Specifically, the anharmonic character of the host-lattice vibrations and local symmetry breaking, which cannot be captured by harmonic phonon calculations alone, have motivated the use of finite-temperature MD simulations.36,72,109,118 However, computational costs associated with doped scenarios of, e.g., Na3PnCh4 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 Na3−xSb1−xWxS4, which was shown to attain excellent ion conductivities.115,116 The mechanism involves W substituting a Sb site, forming [WS4]2− polyanions that exhibit different charges than the original [SbS4]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 Na3SbS4 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 Fig. S10, ESI†). This results in precise predictions of g(r) among all elements in the anion host lattice and the Na+ ions (see Fig. S11, ESI†). 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 vibrations (see Fig. S12, ESI†). The latter suggests the presence of anion–cation coupling, in line with previous studies.28,36,109 Notably, 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 Fig. S13–S15, ESI†), as previously reported.110
Next, we investigate the accuracy of MLFFs for the tungsten-doped Na2.94Sb0.94W0.06S4. A new set of force fields is trained starting from the MLFF of the pristine Na3SbS4. The BEEF values during the training procedure largely remain below 0.02 eV Å−1 after around 1000 steps (see Fig. S16, ESI†). As expected, the error of the MLFF is higher than the one for the pristine system due to increased variations in the local coordination of S atoms in the W-substituted polyanions. Explicit calculation of errors from MLFF to DFT forces were 14.0 meV Å−1, 41.4 meV Å−1, 24.3 meV Å−1 and 52.3 meV Å−1 for Na, Sb, S, and W, respectively (see Fig. S17, ESI†). Comparing these results to the above findings, these forces are in the range of what we found for AgI and below the ones reported for Ge and P in LGPS. Consequently, the MLMD simulations accurately describe g(r) within the [SbS4]3− polyanions as well as between Na+ and S2− ions when compared to AIMD data (see Fig. 8a). We observe only minor deviations in the W–S g(r), where MLMD results show slightly lower peak intensities than the AIMD ones (see also Fig. S18, ESI†). The W–S g(r) exhibits a pronounced peak at 2.21 Å, confirming that bond lengths are shortened compared to Sb–S of 2.36 Å in the pristine case because of the more positive W6+ ions. g(r) for Na–Na remains largely unaffected by doping and still shows a first-order peak at 3.61 Å with broadened long-range distributions.
An important question is whether MLMD simulations accurately describe the vibrational properties of W-doped Na2.94Sb0.94W0.06S4. Indeed, we observe good agreement between MLMD and AIMD data for the overall shape of the VDOS and individual element-wise 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 (Fig. S19, ESI†). 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. (3)), because it was recently identified as a good indicator for changes in the overall lattice stiffness.72 Here, 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 Na3SbS4 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 DNa+ = (1.13 ± 0.69) × 10−6 cm2 s−1. The MSD computed from the single AIMD trajectory lies within the standard deviation of MLMD values. Van Hove correlations for Na2.94Sb0.94W0.06S4 exhibit a behavior that significantly deviates from AgI and LGPS: Gs 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 S20, ESI†). Notably, the region between the two strongly correlated Gs features is only weakly populated, which points to the migration of Na+ ions via thermally activated hopping. Concerted migration can be ruled out from Gd (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 SbS4 and WS4 tetrahedra in the host lattice of Na2.94Sb0.94W0.06S4. Distributions of θ angles within SbS4 and WS4 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 SbS4 tetrahedra exhibit an asymmetric distribution with an extended tail toward larger angle values whereas θ of WS4 tetrahedra exhibit an essentially symmetric characteristic. This indicates the existence of anharmonicity in the rotational motions of SbS4 tetrahedra. This is also visible in the angular distribution, ρ(ω), of SbS4 tetrahedra, which shows small deviations around −0.5 and 0.5 ° fs−1 from the otherwise harmonic behavior. Notably, ρ(ω) of WS4 tetrahedra shows a more harmonic distribution. The overall small deviations between AIMD and MLMD are found to be more pronounced for WS4, which likely is a result of the small sampling size of WS4.
We now focus on correlations between SbS4 and WS4 tetrahedra and neighboring Na+ ions. We obtain the 2D probability density ρ(rNa–S, ϕ), by defining the distance, rNa–S, between neighboring Na+ ions and S atoms, and the dihedral angle ϕ, see Fig. 10 for MLMD and Fig. S21, ESI,† for AIMD results. The largest correlation between SbS4 tetrahedra and Na ions exists in a small range of 80° ≤ ϕ ≤ 100° with rNa–S remaining below ≤4 Å. Minor intensities are found in two regions: one in a range of 100° ≤ ϕ ≤ 120° and rNa–S of ≤4 Å, and the other in 60° ≤ ϕ ≤ 80° and 5 Å ≤ rNa–S ≤ 6 Å. The noticeable intensities at lower rNa–S values showing a broadened ϕ distribution between 80° and 120° indicate that Na+ ions experience a significant degree of rotational freedom around SbS4 tetrahedra. In regard to WS4 tetrahedra, the intensity in the range 60° ≤ ϕ ≤ 80° and 5 Å ≤ rNa–S ≤ 6 Å is found to be more pronounced than in case of SbS4 ones. This observation suggests a facilitated ion migration when Na+ cations are passing by WS4 tetrahedra, which can be explained by W–S bond lengths being shorter than the Sb–S ones and the fact that WS4 tetrahedra occupy less space. Consequently, Na+ ions are 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 as visible in the Gs for Na2.94Sb0.94W0.06S4 (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 Na2.94Sb0.94W0.06S4 are correlated with cation diffusion and suggest that ion conduction can be tuned through adjusting the dynamics of the host lattice.
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,120 Such concerns have led to the development of novel MLFF methods that correct for nonlocal effects including long-range electrostatic interactions.121–124 In the case of SSICs, the short-range nature of interactions and the large dynamical distortions exhibited by the ions likely aid 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 promising approaches towards high-throughput Raman calculations and characterization have been made recently,125–129 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 entirety of vibrational anharmonicity in such materials,130–132 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 finite-temperature Raman spectra viable and inspire novel data-driven approaches, combining ML-models for Raman activities133–139 with MLMD to launch a flourishing synergy between experiments and computation.
Moreover, we envision that computationally efficient and precise MLMD simulations may enrich the physical information that is contained in datasets for a tailored design of novel materials.140–142 Here, the majority of existing approaches for data-driven screening of suitable material candidates mainly relies on properties that are derived in the ground state.14,143–146 Incorporating vibrational properties as well as information about the ion-conduction characteristics predicted from MLMD simulations may improve the selectivity of existing screening approaches and potentially deepen our understanding of the structure–property relations of SSICs.
Footnote |
† Electronic supplementary information (ESI) available: Additional details of the MLFF theory; detailed numerical parameters and settings for AIMD and MLMD simulations; additional figures on the analysis of force error, as well as of structural and dynamical SSIC properties. See DOI: https://doi.org/10.1039/d4ta00452c |
This journal is © The Royal Society of Chemistry 2024 |