Machine-learning-assisted discovery of lattice dynamics signatures of sodium superionic conductors

Abstract

Sodium superionic conductors are key to the development of all-solid-state sodium batteries. Discovery of new superionic conductors has traditionally relied on insights from material defect chemistry and the transition/hopping theory, while the role of lattice vibrations, i.e., phonons, remains underexplored. We identify key lattice dynamics signatures that govern ionic conductivity by analyzing the phonon mean squared displacement (MSD) of Na⁺ ions. By high-throughput screening of a dataset of 3903 Na-containing structures, we establish a strong positive correlation between phonon MSD and diffusion coefficients, providing a quantitative correlation between lattice dynamics and ion transport. To accelerate this discovery, we incorporate machine learning (ML) into our screening workflow, using phonon-derived descriptors to rapidly predict ionic transport properties across a broad structural space. Our findings reveal that low acoustic cutoff phonon frequencies, low center vibrational density of states of Na⁺ ions, slightly higher than the acoustic cutoff frequencies, and enhanced low-frequency vibrational coupling between Na⁺ ions and the host sublattice promote superionic conductivity. Phonon mode analysis further demonstrates that only a small subset of low-frequency acoustic and optic modes contribute dominantly to large phonon MSDs and Na⁺ ion migration, while higher-energy modes contribute negligibly. These insights enable the integration of lattice dynamics descriptors, phonon MSD, Na⁺ VDOS center, acoustic cutoff frequency, and low-frequency phonon coupling into machine learning frameworks, accelerating the discovery and rational design of high-performance sodium superionic conductors.

---

New concepts

Sodium superionic conductors are vital for all-solid-state sodium batteries, yet the role of lattice vibrations in ion transport remains underexplored. By screening 3903 Na-containing structures, we reveal a strong correlation between Na⁺ phonon mean-squared displacement (MSD) and diffusion coefficients, establishing a quantitative link between lattice dynamics and ionic conductivity. Machine learning with phonon-derived descriptors enables rapid prediction across broad structural spaces. We find that low acoustic cutoff frequencies, suppressed Na⁺ vibrational density near the acoustic limit, and enhanced low-frequency coupling with the host lattice promote superionic behavior. Crucially, only a few low-frequency modes dominate Na⁺ migration. These results highlight lattice-dynamic descriptors as powerful tools for machine-learning-driven discovery and design of high-performance superionic conductors.

1 Introduction

Superionic conductors are materials that enable large-scale movement of ions within their lattices, resulting in remarkably high levels of ionic conductivity (σ ≈ 1 mS cm⁻¹ or above at room temperature), which can be observed in liquids, and even in the solid state,^1–3 unlike conventional solids that typically have a low ion mobility. This characteristic renders them extremely valuable for applications as solid-state electrolytes in all-solid-state batteries (ASSBs),⁴ thermoelectric materials,⁵ solid oxide fuel cells,^2,6 and energy storage.⁷ Due to their low manufacturing cost, greater abundance, and comparable chemical and electrochemical properties, Na-ion batteries are a promising alternative to Li-ion batteries.⁸ As a result, sodium superionic conductors (NASICONs) are gaining significant interest as solid-state electrolytes for all-solid-state sodium batteries,³ which are safer when compared to the ignition-prone nature of liquid electrolytes.^9,10 However, there are only a limited number of materials exhibiting exceptionally high ionic conductivity compared to liquid electrolytes,¹¹ which could significantly hinder the big prospects and promising applications of all-solid-state sodium-ion batteries.¹² Therefore, the discovery and design of novel sodium superionic conductors remain an ongoing challenge.³

To address it, computational design methods have been largely employed as a powerful tool for predicting and optimizing superionic conductors.^3,13,14 Techniques such as density functional theory (DFT), ab initio molecular dynamics (AIMD) simulations, and, more recently, machine learning algorithms allow for screening new materials with high ionic conductivity.^15–20 In identifying promising Li SICs, first principles calculations and high-throughput computational screening have already been widely applied due to the increasing computational power and their ability to accurately predict both structural and dynamic properties.² Fujimura et al. identified promising Li SICs with a first principles approach, such as Li₄GeO₄ and Li₄SiO₄, exhibiting high ionic conductivity at elevated temperatures.²¹ Additionally, LiAlSO was predicted to exhibit fast ionic conduction by Wang et al.²² using first principles calculations. Similar approaches are now being adopted for sodium-based systems. An intuitive and straightforward way is to substitute lithium with sodium in Li-SICs, considering that both Li and Na belong to the same alkali metal group in the periodic table and thus their chemical compounds could have similar properties. This approach has been explored in many studies with some success.^3,23,24 While the computational techniques developed for Li SICs lay an invaluable foundation for designing and discovering new sodium SICs, differences in the ion size, mobility, and chemical interactions, all of which impact material structures and ionic transport mechanisms, prevent these methods from being fully applicable for exploring new NASICONs.³ In fact, it has been shown that the knowledge derived from known Li-SICs cannot be simply utilized or transferred to design and discover new Na-SICs,³ despite the chemical similarities between Li⁺ and Na⁺. Moreover, although DFT and AIMD are highly accurate in studying ion transport in NASICONs,^25,26 they are computationally expensive and time-consuming,² which limits their wide deployment in screening of large-scale unknown materials. This limitation opens significant opportunities for leveraging machine learning (ML) to accelerate the discovery of NASICONs.

ML has significantly transformed materials discovery in recent years, offering powerful new tools to predict, design, and discover novel materials with tailored properties faster and more efficiently than traditional approaches. Although ML methods can efficiently screen materials,^17,18,27,28 their effectiveness depends on the availability of large, accurate and diverse datasets to train the ML models, which allow to learn patterns and subsequently make reliable predictions. However, for ion transport in NASICONs, the scarcity of high-quality, diverse training data that cover broad chemical compositions and material symmetries currently poses a grand challenge,²⁹ limiting ML's ability to generalize and accurately predict novel superionic conductors. For instance, the state-of-the-art universal machine learning potentials (MLPs) cannot be straightforwardly deployed for quickly screening fast ionic transport materials. This is understandable considering that the majority of the training data used for training such MLPs are near equilibrium, while ion migration in superionic conductors usually involves out-of-equilibrium movement and thus cannot be easily captured very well by the existing models. Therefore, searching for more efficient material descriptors that have a strong correlation with superionic transport behavior is urgent for the disruptive development and design of novel superionic conductors, either through high-throughput DFT calculations or by machine learning model screening.

In principle, superionic conductors are solid materials in which a subset of the guest atoms (Na⁺ cations herein) can flow as though they are in a “liquid” formed by the host material (such as the anions). From a physical point of view, there is a strong correlation between ionic transport properties and the lattice dynamics (phonons) features of the host materials. Phonons are defined as collective excitations in condensed matter that quantize vibration modes in lattices, which are believed to play a critical role in ensuring structural dynamic stability and enabling the cation transport within the anion framework, but unfortunately, such a role has not been fully understood yet. Based on this premise, using a dataset of ∼3903 Na-containing structures, we propose lattice dynamics features as a new route and design principle for sodium superionic conductors. We demonstrate evidence of strong correlation between the diffusion coefficient of Na⁺ ions and lattice dynamics features, such as lattice softness, low migration energy barrier or activation energy, low acoustic cutoff phonon frequencies, and medium to low central frequency in partial density of states (PDOS) of phonons associated with mobile ions, all of which can be quantitatively characterized by the phonon language in physics, in terms of interatomic force constants (IFCs), mean squared displacements (MSDs), and eigenvectors of vibrational modes of Na⁺ ions at elevated temperatures. A comprehensive understanding of the lattice dynamics of NASICONs can give us insights into how thermal vibrations and structural stability influence ionic transport and provide principles for designing new superionic conductors. Adding lattice dynamics into the previous ML framework is expected to significantly enhance ML model capabilities and accelerate the discovery of next-generation sodium-based solid electrolytes.

2 Computational methodologies

In Fig. 1, we present the schematic of the workflow outlining the key processes of this work. First, 3903 Na-containing structures are filtered out from the Open Quantum Materials Database (OQMD)³⁰ and the Inorganic Crystal Structure Database³¹ (ICSD). The structures are then re-optimized by our own DFT calculations. 293 structures are selected for evaluating universal machine learning potentials, and the EquiformerV2 fine-tuned model on the Open Materials Database (OMAT) and Material Project trajectories (MPtraj)³² model is finally chosen for producing all lattice dynamics properties and MD simulations presented in this work. After screening 921 structures with dynamic stability by the EquiformerV2 fine-tuned model on the OMAT and MPtraj model, i.e., free of imaginary frequencies in the Brillouin zone, MD simulations are performed by the same model, and in parallel, their lattice dynamics features are calculated. Finally, the correlation between ion transport properties and lattice dynamics features is extracted and the insights are further analyzed.

Fig. 1 Schematic workflow outlining the key processes in this study. An initial set of 3903 Na-containing structures was extracted from the OQMD and ICSD databases and re-optimized by DFT. From these, 293 structures were selected for evaluating MLPs, with the EquiformerV2 fine-tuned model (OMAT and MPTraj) ultimately chosen for lattice dynamics calculations and MD simulations. 921 dynamically stable structures were subsequently identified and used to compute lattice dynamics features and perform MD simulations. Finally, correlations between Na⁺ ion transport properties and lattice dynamics descriptors were analyzed to extract key insights.

2.1 Initial structure screening

We screened the OQMD and ICSD and filtered out the sodium-containing structures with a non-zero energy band gap and negative formation energy. Using these criteria, we obtained 3903 sodium-containing inorganic crystal structures. Since the electronic structure and total energy predictions are important for ensuring physically meaningful structural parameters and stable ground-state configurations, density functional theory (DFT) was chosen for structural optimization because of its high accuracy and predictive power. Therefore, we re-optimized all 3903 structures by DFT calculations with our own computational parameters, which are generally stricter than those of the OQMD and ICSD. The computational details of the DFT calculations can be found elsewhere.^33–37 Furthermore, to ensure that these structures are dynamically stable, we used the EquiformerV2 fine-tuned on the OMAT and MPtraj MLP model to calculate the IFCs of the reoptimized structures (see details below on how the MLP is selected) and then screened the phonon dispersions. Upon screening, we found 1363 structures to be dynamically stable with positive dispersions, i.e., absent from negative frequencies in the Brillouin zone. The statistical analysis of the elemental distributions of the 1363 Na-containing structures across the periodic table after this screening is presented in Fig. S1 in the SI. The space group number comparison before and after screening for positive dispersions is shown in Fig. S2. Our data are represented by a total of 61 elements spanning the periodic table and all 230 space groups, classified into triclinic, monoclinic, orthorhombic, tetragonal, trigonal, hexagonal, and cubic systems (Fig. S1). We found that monoclinic, orthorhombic, tetragonal, and trigonal space groups constitute the majority of our dataset after the initial screening of positive dispersions.

2.2 Evaluating machine learning potentials

The accuracy of an MLP directly determines the quality of an MD simulation that relies on it. Therefore, a prerequisite step for high-throughput MD simulations and subsequent lattice dynamics analysis is to evaluate the performance of MLP models and then select the best model in terms of accurately predicting atomic forces and potential wells in complex energy landscapes. In this study, we systematically evaluated the performance of various ML models in predicting atomic forces on a dataset comprising 293 Na-containing structures taken from the OQMD. The evaluated MLP models include the EquiformerV2 pre-trained model,³⁸ the EquiformerV2 fine-tuned model trained on OMAT and MPtraj,³² the EquiformerV2 fine-tuned model trained solely on OMAT,³² the MatterSim³⁹ pre-trained model, the MACE⁴⁰ pre-trained model, and the CHGNET⁴¹ pre-trained model. After structure re-optimization, supercells were created by expanding the optimized primitive cells with a different supercell size suitable for the structure. The supercell size is generally determined by fulfilling the following requirements: (1) the lattice parameter of the supercell in different crystallographic directions is similar and (2) the total number of atoms in supercells must be at least 80; indeed, most of the finalized supercells contain atoms with total numbers in the range of 80 to 300. The atoms in the supercells are then displaced randomly in different directions by a constant displacement of 0.03 Å using the PHONOPY⁴² package. By doing so, 30 different displaced supercells are generated for each structure, and the atomic forces in the supercells are calculated using high-accuracy self-consistent DFT. All 293 structures have 30 displaced supercells each, making a total of 8790 datasets, which are used to evaluate the different MLP models. The predicted atomic forces obtained from these MLP models were benchmarked against forces computed directly by DFT. Root mean square error (RMSE) values between the MLP predicted forces and the DFT calculated forces were computed to quantitatively assess the performance of each MLP model. Our results indicated that the EquiformerV2 fine-tuned model trained jointly on OMAT and MPtraj exhibited one of the lowest RMSE values (57.7 meV Å⁻¹) when compared to DFT forces, though not the absolute lowest among the total of 6 MLPs. However, as discussed below, we chose the EquiformerV2 fine-tuned model in our subsequent analysis based on additional factors such as its superior performance for predicting the interatomic force constants (IFCs) and phonon stability.

While the quality of atomic force prediction is very important for an MD simulation, the derivatives of interatomic forces, in particular the 2nd order, are crucial for our insight analysis, such as the focus of MSD in this study. To this end, we also compared the 2nd and 3rd order interatomic force constants (IFCs) derived from each ML model against their corresponding IFCs calculated by DFT to further validate and strengthen our model selection. When it comes to accurately capturing anharmonic effects and lattice dynamics, the IFCs derived from DFT computations serve as an essential benchmark because of their higher predictive power and physical rigor. Consistently, we find that the EquiformerV2 fine-tuned model on OMAT and MPTraj³² demonstrated the lowest RMSE values (367.8 meV Å⁻² and 8163 meV Å⁻³) when compared against DFT for both 2nd and 3rd order terms (Fig. S3 and S4). Furthermore, phonon dispersion analysis revealed that all 293 Na-containing structures exhibited positive dispersions by DFT. Upon evaluating the MLP-derived dispersions, the EquiformerV2 fine-tuned model on OMAT and MPtraj yielded positive dispersion for 258 structures, i.e., a success rate of 88%, surpassing all other ML models in terms of dynamic stability prediction accuracy. In light of the EquiformerV2 fine-tuned model's exceptional capacity to precisely reproduce atomic forces and IFCs as benchmarked against DFT calculations, as well as providing the most reliable phonon stability predictions, we finally selected this model for further MD simulations and lattice dynamics analysis throughout this study. The demonstrated accuracy and robustness of the EquiformerV2 fine-tuned ML model make it particularly suitable for investigating lattice dynamics and related transport properties in complex materials. Additionally, after screening all promising Na-superionic conductors, we re-examined the performance of the top 3 MLP models, i.e., EquiformerV2 fine-tuned on OMAT and MPtraj, EquiformerV2 fine-tuned on OMAT only, and the MatterSim pre-trained model, by randomly sampling 200 atomic configurations from each AIMD run of 31 Na-superionic conductors and compared the force evaluation performance in Fig. S5 in the SI. The results indicate that the EquiformerV2 model fine-tuned on MPtraj and OMAT achieved the best performance with an RMSE of 11.6 meV Å⁻¹ as compared to AIMD forces. This further underscores our decision to employ this model for the remainder of this work.

We also would like to point out that the field of universal machine learning potential, in particular for materials science, has developed very fast in recent years. Lots of new MLPs come out just every few months. When we initiated this work, the EquiformerV2 model was in the very top position of the MLP list available in the Matbench list, which is the major reason why we include it in our ML model test. We acknowledge that currently Matlantis⁴³ is probably one of the most accurate MLPs in materials science. However, as the MLP algorithms and methodologies iterate very fast, and the high-quality DFT data used for training such models are also growing very fast, one cannot always catch up with the latest model trained on the newest data and deploy it in their study. Systematically evaluating very large numbers of MLPs including the newest development of MLPs for ionic transport purposes would be a topic of future studies.⁴⁴

Regarding ML model's generalization problem, for our diffusion coefficient calculations by MD simulations, the quality of atomic forces is extremely important, instead of energy prediction, as the atomic forces directly determine how the atoms/ions are moving (trajectories) in the lattice according to Newton's second law of motion, from which the diffusion coefficient is determined by the Einstein relation. Moreover, our idea of quantifying and correlating superionic behavior with lattice dynamics requires accurate prediction of the 2nd order force constants, which are directly related to the quality of atomic forces, since physically speaking the 2nd order force constants are the 2nd derivative of potential energy with respect to atomic displacement and therefore can be regarded as the 1st derivative of atomic forces with respect to atomic displacement as well.

The OMAT24 dataset, in which the EquiformerV2 is fine-tuned, contains over 110 million diverse, far-from-equilibrium DFT structures, thus providing a wide range of chemical space, while the MPtraj dataset, obtained from the materials project relaxation trajectories, improves the models’ accuracy on equilibrium structures.³² This combination enables the EquiformerV2 fine-tuned model to achieve accuracy in atomic forces, positioning it well-suited for this study. Therefore, we believe that the EquiformerV2 model fine-tuned on the MPtraj and OMAT dataset provides great generalization ability for predicting decent quality of atomic forces in various atomic environments, such as those Na-structures studied herein. This is reflected by the low RMSE values of predicted forces as compared with AIMD calculations, as evidenced in Fig. S5.

2.3 MD simulation with EquiformerV2 fine-tuned machine learning potential

For each primitive cell of the 1363 structures, we created a supercell with the same size as the previous dynamic stability screening using the PHONOPY⁴² package. The supercells were then used for performing classical MD simulations with the EquiformerV2 fine-tuned MLP.³² The MD simulations were run with a constant volume and constant temperature (NVT) ensemble at T = 800 K with a Langevin thermostat. Notice that elevated temperatures will facilitate the movement of Na-ions within the supercells while preserving the original crystalline structures, though excessively high temperatures might melt or destroy the structures. The MD simulation for each structure consists of two stages: (1) the first 10 000 steps are for structure relaxation and (2) the subsequent 40 000 steps are for a production run with trajectory output every 4 MD steps for post-processing of the diffusion coefficient. A timestep of 1.5 fs is used throughout.

2.4 Post-processing MD data

Mean-squared displacement (MSD) data are retrieved from the trajectories of the EquiformerV2 fine-tuned MD run. The MSDs of each constituent element in a specific structure are computed based on the time-dependent atomic positions in all directions (x, y, z), which can be expressed aswhere N is the number of all atoms belonging to that element in the system, r_i(t) is the position of atom i at time t, and r_i(0) is the initial position of atom i. Calculating the MSDs provides insights into the diffusive behavior of the atoms in the simulated system. The diffusion coefficients (DCs) of each constituent element in the structure were calculated bywhere n is the dimensionality of the system (herein, n = 3 for average DC, while n = 1 for maximum DC). When calculating the diffusion coefficients, the slope of the MSD data is fitted by a linear function with a zero intercept, which is reasonable since there is no diffusion at the beginning (t = 0). It is worth pointing out that we need to filter out melted or destroyed structures during the MD run where the selected MLP may poorly represent those materials, or the materials are not stable at the temperature T = 800 K. To this end, the diffusion coefficients obtained were used in filtering out structures in which the non-Na atoms exhibit high diffusion coefficients on the order of 10⁻² Ų ps⁻¹ or above in all three spatial directions. In addition, we implemented a method of calculating the average of the MSD data for the last 20% period of the MD production run. Using this additional criterion, structures containing non-Na atoms with a high average MSD greater than 1 Ų in any of the three spatial directions are duly filtered out. This method ensures that all the unstable structures are filtered out, bringing our total structures from 1363 down to 921, from which all the reported results herein are calculated. Our approach is an effective means to identify and eliminate unstable structures with high ionic mobility of non-Na atoms across all spatial directions (x, y, and z).

2.5 Phonon calculations using the EquiformerV2 fine-tuned model

For the phonon calculations, following structure re-optimization by DFT, we used the PHONOPY⁴² package to generate randomly displaced supercells, each with an atomic displacement magnitude of approximately 0.03 Å. This approach is advantageous as it systematically probes the local energy landscape, allowing for accurate determination of interatomic force constants (IFCs), which characterize the local potential wells and energy barriers between neighboring potential wells. Subsequently, we evaluated the atomic forces on the displaced supercells using the selected EquiformerV2 fine-tuned MLP model. These forces serve as the input to fit both the harmonic (2nd order) and anharmonic (3rd order) IFCs by the compressive sensing lattice dynamics (CSLD) method.^45–47 The forces and displacement configurations are used to construct the sensing matrix A and the force vector F. The IFCs are determined by solving the following optimization problem:where A is the sensing matrix constructed from atomic displacements. The second term is the standard Euclidean [small script l]₂ norm, which measures the force-fitting error over the training data. The first term, the [small script l]₁ norm, promotes sparsity by driving a small number of nonzero IFCs. μ is the weight adjustment term for [small script l]₁ and [small script l]₂ terms, controlling the trade-off between sparsity and fitting accuracy. By promoting sparsity and accuracy, this approach efficiently captures relevant lattice interactions, making CSLD highly suitable for modelling strongly anharmonic systems with minimal computational cost. We used the Pheasy package, which implements the CLSD algorithm to do the IFC fitting.⁴⁸ The second-order force constant is given by⁴²where V is the potential energy, r_α(i) is the displacement of the ith atom in the direction α and r_β(j) is the displacement of the jth atom in direction β. The IFCs determine the bond strength in the lattice and have a direct impact on the thermal properties and ionic transport of the material, as we will see later.

We then used 2nd-order IFCs to calculate the MSD using the PHONOPY⁴² package. The MSD calculated based on IFCs is crucial for understanding thermal vibrations and atomic diffusions, which provides deep insight into how the atoms move away from their equilibrium position over a specific time interval due to thermal effects. The MSD is calculated as follows:⁴²

where j and l are the labels for the jth atomic position in the lth unit cell, t is the time, α is the Cartesian axis, ℏ is the reduced Planck constant, N is the number of unit cells, m is the atomic mass, q is the wave vector, v is the phonon mode index, e is the polarization vector of the atom (j,l) and the phonon band q, and n_v(q,T) is the phonon population given bywhere k_B is the Boltzmann constant and T is the absolute temperature.

To examine the contribution of the various phonon modes to Na⁺ ion displacement, we computed the MSD band by band by decomposing the general MSD equation [eqn (5)] for each phonon mode at every q-point sampled in the Brillouin zone. We used the mode frequency and the amplitudes of the eigenvectors projected onto the Na⁺ ions with the harmonic approximation as implemented by PHONOPY.⁴² For each phonon mode, the MSD was computed by summing the components of the squared eigenvectors over the three Cartesian directions and weighting with the Bose–Einstein occupation factor and inverse frequency. This decomposition allows one to identify phonon mode level contributions to the large displacement of Na⁺ ions and demonstrates the lattice dynamics origin of ionic mobility, such as some low-frequency phonon modes that are accountable for the large vibrational amplitudes of Na⁺ ions.

We calculated the center of the vibrational density of states (VDOS) for Na⁺ ions from the MD trajectory through fast Fourier transformation of the autocorrelation function of velocities. The VDOS center, which represents the first moment of the Na-specific vibrational spectrum, is defined as the average frequency weighted by the VDOS^49,50

where ω_i is the ith vibrational frequency, and g(ω_i) is the corresponding VDOS value at that frequency. A lower VDOS center suggests that the lattice is soft, with low-frequency vibration, which facilitates ion mobility. In comparison, a higher VDOS center suggests a stiffer, more rigid lattice, which potentially restricts ion mobility.

Furthermore, we computed the phonon participation ratio (P_λ) for each mode to assess the degree of mode localization. The participation ratio for a specific phonon mode λ is expressed as⁵¹

where N is the number of atoms in the primitive cell, e_ia,λ is the eigenvector along the α direction (x, y, and z) of the ith atom in the primitive cell for each vibrational mode at q and band index v and * is the complex conjugate operation. A high participation ratio value (close to 1) shows that the phonon mode is delocalized over many atoms, meaning almost all atoms in the primitive cell participate in the vibration of the mode, while a low participation ratio value (close to 0) denotes phonon localization where only a small portion of atoms and even a single atom participate in the vibration of the mode.⁵¹

To quantify the correlated vibration between mobile Na⁺ ions and the surrounding host sublattice, we define a phonon mode resolved Na-host coupling factor for each phonon mode (ω, ν):

where e⃑_Na⁺,α and e⃑_host,α are the phonon eigenvector components along the Cartesian direction (α = x, y, z) of mobile Na⁺ ions and the host sublattice, respectively. Coupling factors close to +1 (−1) indicate strongly correlated in-phase (out-of-phase) motion between Na⁺ ions and host sublattices, while coupling factors close to 0 indicate uncorrelated motion.

2.6 AIMD simulations and phonon calculations by DFT

While the EquiformerV2 fine-tuned model jointly trained on the OMAT and MPtraj datasets is efficient in providing fast evaluation of interatomic forces, its accuracy is limited by the quality and range of data it was trained on. Therefore, some uncertainty may arise when configurations fall outside its training dataset, which will affect both the phonon bands through the random displacement technique and the diffusion coefficient results through MD trajectories. To ensure the reliability of our results, in particular the trend of the dependence of Na⁺ ion diffusion on the phonon features, we employed AIMD simulations for selected structures to validate the diffusion coefficient results from the EquiformerV2 fine-tuned model. Γ-Point is used for electrons in all AIMD simulations with a cutoff energy of 520 eV and an energy convergence criterion of 10⁻⁶ eV. The total simulation time ranges from 30 to 90 ps or even longer with a timestep of 1.5 fs. All AIMD simulations were run with the NVT ensemble at 800 K and with the Nose–Hoover thermostat. The time-dependent MSD curves for four selected structures with predicted high DC by AIMD are provided in Fig. S6 in the SI, with fitted diffusion coefficients labeled. In addition, we directly compared the time-dependent MSD curves obtained from AIMD and the EquiformerV2 fine-tuned MLP for the four selected structures in Fig. S7. Our results of the MSD curves indicate a good agreement between our preferred MLP model and AIMD. Furthermore, the phonon band structures of the four representative materials are further compared between the EquiformerV2 fine-tuned model and DFT. As presented in Fig. S8, the phonon dispersions by the two methods are essentially the same, again confirming the reliability of the EquiformerV2 fine-tuned model.

3 Results and discussion

3.1 MSD results by supercell calculation

Although our primary analysis is based on the MSD derived from the primitive cell for our calculations due to computational efficiency and ease of screening numerous structures, it is essential to note the inherent limitations of this method. Primitive cells capture the dynamics based solely on the globally minimized energy configurations at 0 K, potentially underestimating realistic ionic transport behaviors at elevated temperatures. To illustrate the significance of employing the supercells, Fig. S9a and b show the results of the average and maximum MSD, respectively, of a representative structure, Na₄TeSe (OQMD ID 1554824), calculated by the supercell sampling at different random snapshots from the AIMD trajectory. First, it is worth pointing out that supercells are necessary for calculating the realistic MSDs of the superionic conductors with moving ions inside. This is because primitive cell calculations only use a single atomic configuration (most likely global minima of the structure obtained by DFT at 0 K) and thus only yield a specific MSD value for a given structure. Such a single potential energy minimized scenario cannot fully represent the realistic dynamically changing configurations of ionic systems at elevated temperatures, in particular when the partial ions become movable and migrate far away from the original global minima positions. As a result, primitive cell calculation often underestimates the MSD values since the globally minimized configurations usually correspond to deep potential wells and thus the local bonding strength is much stronger than the “randomized” local minima with a prevalent loosely bonded atomic environment (see details below). In Fig. S9a, the small deviation across the different time snapshots suggests that when the MSD is averaged over all Na⁺ ions in all three directions, the structure exhibits consistent phonon dynamics features governing the ionic movement. This implies that the overall thermal motion of Na⁺ ions is balanced and stable, indicating uniformity in its ionic bulk properties, and the underlying dynamics of the system are stable across the different snapshots. However, this does not mean the lattice vibrations or local atomic environments of Na⁺ ions residing on different sites are the same. In contrast, some Na⁺ ions at some snapshots show distinctly different levels of ionic motions. The large deviation across the different time snapshots shown in Fig. S9b is indicative of exceptionally large thermal displacements experienced by some individual Na⁺ ions in different spatial directions. These individual Na⁺ ions with high MSD play a significant role in contributing to the overall ionic conductivity of the structure. It is also important to note that there are usually only a small portion of Na⁺ ions possessing high MSD at a certain snapshot, as shown in Fig. S9c. Among all 72 Na⁺ ions in this structure, only 5 ions have high MSD well above the average level. This can be explained in terms of those ions having a more favorable environment, for example having fewer neighbors or loosely bonded with anion neighbors, thereby exhibiting much higher MSD values. The dynamic change of all MSD results in Fig. S9 provides strong evidence of the importance of using supercells with randomized configurations of moving ions extracted from MD trajectories, as this approach captures realistic dynamics in superionic conductors. However, due to computational constraints, our present study employs the use of the primitive cell, recognizing that while the primitive cell might underestimate the true MSD value compared to the supercell, using primitive cells still provides valuable insights for efficiently identifying correlations between the MSD and diffusion coefficients across a large set of structures.

3.2 Lattice dynamics features for sodium ion transport

Fig. 2(a) and (b) demonstrates the average and maximum diffusion coefficients of the 921 structures against their MSD, respectively, i.e., our proposed lattice dynamics features. Note that all diffusion coefficient values are obtained by post-processing MD trajectories of supercells by the EquiformerV2 fine-tuned model, unless explicitly specified from AIMD runs otherwise. The MSD of the structures in Fig. 2 is obtained from the primitive cell calculation. We observe a positive correlation between the MSD and the diffusion coefficients. This trend implies that structures with high MSD, which is indicative of large random displacements, tend also to exhibit high diffusion coefficients. This is in line with the Einstein relation given by D ∝ 〈u²〉 and with well-established theories that high vibrational amplitudes are frequently associated with enhanced ion mobility in superionic conductors.^50,52 However, it can be observed that some structures deviate from this trend, with relatively high DC and low MSD and vice versa. This can be attributed to several factors. Certain materials show high MSD due to localized rattling of ions within small cages instead of the actual long-range migration. This is common in soft or disordered frameworks, such as argyrodites and thiophosphates, where local vibrational motion predominates but makes a negligible contribution to new diffusion.⁴⁹ Long-range ion transport can be impeded by large energy barriers between sites or insufficient connections between diffusion channels, even in cases where ions are highly mobile and show active vibrations. Despite considerable local motion, this results in low or modest diffusion coefficients. To further validate our results using the MLP method, we selected a subset of our data with 36 representative structures comprising both low and high diffusion coefficients and then ran AIMD to obtain diffusion coefficients with MSD values calculated by DFT (red color in Fig. 2) based on 2nd order IFCs using the same primitive cell. As established earlier, DFT accurately captures the IFCs and phonon spectra, which directly influence thermal displacement amplitudes. Comparison between DFT and MLP shows strong consistency, thereby confirming the accuracy of the EquiformerV2³² model in capturing vibrational features relevant to ionic transport. While EquiformerV2 fine-tuned model equipped MD is computationally efficient for calculating the diffusion coefficients of considerably large number of supercells with diverse material symmetry and compositions (a typical 200-atom supercell MD run of total 50 000 steps only takes about 3 hours on one GPU node with the EquiformerV2 fine-tuned model), the accuracy of the diffusion coefficients can be limited by uncertainties in the predicted atomic forces. In contrast, AIMD offers greater accuracy in capturing the true atomic forces, leading to more accurate diffusion coefficients.

Fig. 2 (a) Average and (b) maximum diffusion coefficients of Na⁺ ions against the corresponding mean squared displacement (MSD) of 921 Na-structures with MSD calculated using primitive cells and EquiformerV2 fine-tuned model on OMAT and MPTraj (blue circles), 36 Na-structures with diffusion coefficients by AIMD and MSD validated by DFT (red squares), and 7 previously reported structures^54,69–71 (magenta stars). All results are calculated at 800 K. The formulas of four representative structures are indicated by arrows. (Inset) Schematic of the comparison of different migration energy barriers between small and large MSDs. The orange arrow indicates the softer lattice when the MSD becomes larger.

In Fig. 2(a) and (b), we see a clear trend that both the average and maximum diffusion coefficients increase with their corresponding MSD. This strong correlation or lattice dynamics signature of ion diffusivity can be explained by the underlying link between both the Na⁺ ions' thermal motion behavior and migration with their local potential well or energy landscape. Statistically speaking, most of the time, the Na⁺ ions stay in the local potential well and simply do random thermal motion there. However, they have some chances to escape the local potential well at some point and then migrate or jump to the neighboring local potential well. Such probability is governed by two major factors: (1) how far the Na⁺ ions vibrate from their equilibrium positions with given thermal energy, which roughly scales with , where T is the system temperature and (2) how high the energy barrier between two neighboring potential wells is. Actually, both can be reflected by the MSD. On one hand, a higher MSD suggests that the potential wells in which the Na⁺ ions reside are very flat, such that the Na⁺ ions can move far away from equilibrium positions at a given temperature and thus they have higher chances to migrate to the neighboring sites. On the other hand, a higher MSD corresponds to a softer lattice^52–54 or loose bonding, so the energy barrier between two neighboring sites with a certain distance would be much lower (see schematic in the inset of Fig. 2(b) for comparison between high and low MSD cases). This phonon feature of high MSD can therefore be used for fast screening structures exhibiting high Na⁺ ion mobility for sodium all-solid-state battery applications.

Following our analysis of the correlation between the DC and MSD, we further investigated the activation energy barrier of Na⁺ migration for the four representative structures with high DC predicted by our MLP model and validated by AIMD simulations as exhibiting superionic behavior for a temperature range of 300 to 800 K. The activation energy is obtained from the temperature-dependent diffusion coefficient from MD simulations with the EquiformerV2 fine-tuned MLP by fitting the Arrhenius relation , where D is the diffusion coefficient, D₀ is the pre-exponential factor, E_a is the activation energy, k_B is the Boltzmann constant, and T is the absolute temperature. The E_a values for the representative structures, ranging from 133 to 198 meV and reflecting the depth of the potential wells that govern Na⁺ migration, are presented in Table 1. These results fall within the range of reported Na-superionic conductors, such as γ-Na₃PS₄ (161 ± 13 meV),⁵⁵ c-Na₃PSe₄ (280 meV),⁵⁶ Na₁₀SnP₂S₁₂ (350 meV),¹³ and Na₃SbSe₄ (193 meV),⁵⁷ with two candidates (Na₃YBr₆, OQMD #1750683 and Na₃ScBr₆, ICSD #401335) exhibiting barriers comparable to or even lower than γ-Na₃PS₄. This underscores the effectiveness of the descriptors employed in our screening approach in capturing low-energy barriers as fast ion transport characteristics. These results also underscore the capability of the EquiformerV2 fine-tuned model to accurately predict superionic behavior, thus enabling the efficient identification of promising Na-superionic conductors.

Table 1 Activation energy of the four representative structures calculated using the EquiformerV2 fine-tuned model

Database	Material ID	Formula	Energy barrier (meV)
OQMD	1391859	Na4TeS	198.2
OQMD	1554824	NaVPdS4	192.3
OQMD	1750683	Na3YBr6	151.5
ICSD	401335	Na3ScBr6	133.1

---

3.3 Correlation between other phonon features and diffusion coefficients

We have established the general trend and relationship between diffusion coefficients and MSD, one of the most important phonon features of the lattice. It is noteworthy to explore whether other phonon features correlate with the dynamics of Na⁺ ions, together with the traditionally widely used simple structural descriptors. In Fig. S10a, we show the Pearson correlation between the material descriptors and the maximum and average diffusion coefficients. A correlation of 1 indicates a perfect positive relationship, which means as one variable increases, so does the other, while a correlation of −1 indicates a perfect negative relationship, which indicates that as one variable increases, the other decreases.⁵⁸ The traditional structural descriptors shown in Fig. S10a offer insight into how the atomic arrangement of a material affects ion mobility. These descriptors can directly influence diffusion pathways of Na⁺ ions by their interaction with their local environments. The negative correlation between the average number of nearest neighbors for Na⁺ ions and the maximum and average diffusion coefficients (−0.46 and −0.41, respectively) indicates that a higher coordination number (neighboring atoms surrounding Na⁺) restricts Na⁺ diffusion; thus, fewer neighboring atoms surrounding Na⁺ provide more room for Na⁺ diffusion. The negative correlation with the number density (−0.58 and −0.52) and packing fraction (−0.30 and −0.29) supports this idea, because a higher number density implies a more tightly packed structure, which restricts ionic movement by reducing the volume for Na⁺ to diffuse. To this note, the volume shows a positive correlation of 0.38 with the diffusion coefficient, indicating that a larger volume is beneficial for fast Na⁺ ion diffusion, which is very intuitive. The negative correlation with Pauli electronegativity (−0.41 and −0.36) indicates that an increase in the overall electronegativity of the surrounding atoms creates a stronger electrostatic interaction that traps the Na⁺ ions, thereby limiting the movement of Na⁺ ions and making it difficult to escape its local atomic environment. The diffusion coefficient has a slight positive correlation with the average electron count (0.26 and 0.23), which indicates that materials with a higher electron density allow efficient Na⁺ ion diffusion through favorable interactions with the electronic structure of the material. The average bond length of all atoms within a radius of 3.5 Å shows a moderate correlation (0.43 and 0.39) with the diffusion of Na⁺, which indicates that longer bond lengths within this cut-off radius favor Na⁺ ion mobility. These correlations emphasize the influence of the material's structural properties on the diffusion of Na⁺ in the material.

As evidenced in Fig. 2, the vibrational dynamics of the materials have a significant impact on the diffusion behavior of Na⁺ ions, since the ion mobility is influenced by thermal vibrations and lattice dynamics of both movable Na⁺ ions and immobile anion frameworks. In Fig. S10b, we present the Pearson correlation of other major phonon features with the maximum and average diffusion coefficients of the materials analyzed on our 921 data. The partial phonon density of state (PDOS) of Na⁺ ions describes the vibrational frequency only contributed by the Na⁺ ions in the material and accounts for specific energies or frequency ranges. The total phonon density of states represents the contributions of all ions in the material. The center PDOS feature is the average phonon frequency weighted by the DOS at which the Na⁺ ions predominantly contribute to the vibrational states of the material. The skew and kurtosis describe the shape of the PDOS and total DOS distribution. Skewness is a measure of how asymmetric the PDOS and total DOS are, while kurtosis describes whether the distribution is more peaked or flat when compared to the normal distribution. The width of the PDOS and total DOS is the range of frequencies in which the Na⁺ ions and all ions contribute to the vibrational modes, respectively. A wider (smaller) width of PDOS indicates that the Na⁺ ions participate in a broad (narrow) frequency range. We first noticed the strongest positive correlation (0.87 and above) between previously calculated MSDs and diffusion coefficients, which is consistent with the results presented in Fig. 2. We then found a strong negative correlation between the PDOS center of Na⁺ ions and the maximum and average diffusion coefficients (−0.68 and −0.75, respectively), indicating that the lower center of PDOS of Na⁺ ions corresponds to lower frequency vibrations, which facilitates ionic diffusion. This implies that the Na⁺ ions diffuse within the materials through participating in low-frequency phonon modes with lower energies. Higher frequency modes indicate strong interatomic bonding and a stiffer lattice, which limits Na⁺ ion mobility. In contrast, the correlation between the center of the total DOS and the maximum and average diffusion coefficients (−0.48 and −0.54, respectively) becomes weaker when all ions in the material are considered, which shows the impact of non-diffusing ions on limiting the mobility of Na⁺ ions. At this moment, it is hard to distinguish the positive or negative impact of non-diffusing ions in a specific material on the diffusivity of mobile Na⁺ ions, which deserves further study on the detailed “entanglement” or coupling between the movable ions and the immobile sublattices. The negative correlation with skew (−0.31 and −0.35) indicates that more symmetrical PDOS distribution favors high Na⁺ ion diffusion. Symmetrical PDOS means that Na⁺ ions participate equally in both low and high-frequency phonon modes. The negative correlation with kurtosis (−0.33 and −0.30) indicates that a sharper and more peaked distribution (corresponding to lower kurtosis) of PDOS supports Na⁺ ion diffusion.

The acoustic cutoff phonon frequency, which is defined as the maximum vibrational frequency among the three acoustic phonon branches (two transverse and one longitudinal) obtained from phonon dispersion calculations, also shows a negative correlation (−0.56 and −0.52) with the diffusion coefficient. This correlation demonstrates that a lower acoustic cutoff frequency is associated with high Na⁺ ion mobility, which is consistent with the idea that a softer lattice allows the free movement of ions.^2,59 The detailed inverse correlation is demonstrated in Fig. S11a and b. A lower acoustic cutoff frequency indicates that the lattice is soft and that the Na⁺ ions can vibrate with large amplitude at lower energy. This implies that the soft lattice environment creates favorable conditions in which the Na⁺ ions experience less resistance as they diffuse through the lattice, thus facilitating fast Na⁺ ion diffusion.¹¹ The significance of this is the creation of a rattling effect in the lattice of the material, with the Na⁺ ions having enough room to rattle about their equilibrium positions, since the lattice does not restrict them tightly. Hence, the increased mobility of the Na⁺ ions allows them to easily overcome the energy barrier, providing a conducive environment for fast Na⁺ ion diffusion, thus contributing to higher diffusion coefficients and ionic conductivity. In contrast, materials with a higher acoustic cutoff frequency exhibit a lower diffusion coefficient due to their lattice being stiffer, which restricts the free movement of Na⁺ ions in the materials, thus resulting in a lower diffusion coefficient of the materials.

Fig. 3(a) and (b) show a stronger inverse correlation between the VDOS center of Na⁺ ions and both the average and maximum diffusion coefficients, compared to the correlation observed with acoustic cutoff frequency in Fig. S11a and b. The steeper slope in Fig. 3(a) and (b) indicates that structures with a lower VDOS center of Na⁺ ions exhibit significantly higher ionic mobility. This suggests that when vibrational modes associated with Na⁺ ions are concentrated at the lower frequencies, diffusion is enhanced,⁶⁰ a trend consistent with that observed at lower acoustic cutoff frequencies but more pronounced. The color bar in Fig. 3(a) and (b) is the ratio of the VDOS center for Na⁺ to the acoustic cutoff frequency of that material. Regarding purple symbols, the central phonon frequency of Na⁺ ions is medium to high but is still comparable to or even less than the acoustic cutoff frequency. Although this indicates that the dominant phonon modes the Na⁺ ions participate in are acoustic, their diffusion coefficients are generally not high because of the relatively high acoustic cutoff frequency, as explained in Fig. S11a and b. The green and blue symbols indicate that the VDOS center of Na⁺ ions in those materials is just above the acoustic cutoff frequency. These materials have some chances to achieve high diffusion coefficients, if they happen to fulfill the conditions of both low acoustic cutoff frequency (around 1 to 2 THz) and low central frequency in VDOS of Na⁺ ions (around 2.5 to 4 THz). This finding seems to be consistent with the previous general thought that low-frequency collective vibrations of the lattice play a critical role in fast ionic transport.^49,50 It should be noted that, since the central frequency of Na⁺ ions is already higher than the acoustic cutoff frequency of the material, this means that the majority of phonon modes the Na⁺ ions participate in are not acoustic, but optic. Another important effect of the central frequency of Na⁺ ions higher than the acoustic cutoff frequency is the higher chance for more matched VDOS in the low frequency region which results in more correlated motion (see more details below). For the magenta symbols and beyond in Fig. 3(a) and (b), where the central frequency of Na⁺ ions is far greater than the acoustic cutoff frequency, these materials generally have mild diffusion coefficients, meaning that the Na⁺ ions do not have a high chance to couple with low-frequency phonon modes in the lattice. Overall, Fig. 3 together with Fig. S11 demonstrates that a low acoustic cutoff frequency and a low central frequency of Na⁺ ions, slightly higher than the acoustic cutoff frequency, are desirable for achieving a high diffusion coefficient in the material.

Fig. 3 The correlation of (a) average and (b) maximum diffusion coefficients at 800 K with the center of partial vibrational density of states (VDOS) of Na⁺ ions of the 921 dynamically stable Na-structures. The color bar indicates the relative ratio of the center frequency of the partial vibrational density of states of Na⁺ ions to the corresponding acoustic cutoff frequency of that material. Lower center frequency of Na⁺ ions but slightly higher than the acoustic cutoff frequency promotes superionic transport.

3.4 Deep insight into phonon mode level lattice dynamics and correlated ionic transport mechanisms

In Fig. 4(a)–(l), we show the phonon dispersion (a, d, g, and j), frequency dependent participation ratio (b, e, h, and k), with both colored by the phonon mode resolved MSD contribution to the Na⁺ atoms, and Na⁺ ion transport channel (c, f, i, l) in the material. Four representative structures with high diffusion coefficients confirmed by AIMD simulation and large MSD are presented: (a–c) Na₄TeS OQMD ID 1554824 (average and maximum diffusion coefficients 1.24 × 10⁻¹ and 1.64 × 10⁻¹ Ų ps⁻¹, respectively), (d–f) NaVPdS₄ (OQMD ID 1391859) (average and maximum diffusion coefficients 1.96 and 6.53 × 10⁻¹ Ų ps⁻¹, respectively), (g–i) Na₃YBr₆ (OQMD ID 1750683) (average and maximum diffusion coefficients 1.00 × 10⁻¹ and 1.26 × 10⁻¹ Ų ps⁻¹, respectively), and (j–l) Na₃ScBr₆ (ICSD ID 401335) (average and maximum diffusion coefficients 8.03 × 10⁻² and 7.18 × 10⁻² Ų ps⁻¹, respectively). All phonon dispersions and participation ratios shown in Fig. 4 were computed by DFT in order to accurately explore the microscopic lattice dynamics features underlying ionic mobility. We analyze these structures using phonon dispersion relations and their corresponding phonon participation ratio, with the decomposed MSD of Na⁺ ion projected as a color map across each band. As observed in the phonon dispersion (left panel in Fig. 4), the low-frequency phonon mode, typically below ∼1.5 THz, contributes dominantly to the MSD of Na⁺ ions as indicated by the red color. Those phonon modes belong to acoustic and low-lying low-energy optical branches. The dominance of the acoustic and low-lying optical phonon mode highlights their role in facilitating the mobility of Na⁺ ions in the lattice. These low-frequency vibrations with large MSD contribution correspond to lattice softening, in which the atoms experience shallow potential wells (see schematic in the inset of Fig. 2) and undergo large thermal displacements from their equilibrium position and consequently induce a lower migration energy barrier and easier hopping of Na⁺ ions,⁶¹ thus resulting in the high ionic conductivity of these materials. In contrast, the majority of phonon modes contribute negligibly to the MSDs of Na⁺ ions, meaning that the Na⁺ ions participating in those phonon modes just do random thermal motion and hardly migrate. This behavior is in line with previous studies in which low phonon frequencies have been associated with superionic conductivity.^{53,54,61–63} The corresponding participation ratio plots (middle panel in Fig. 4) further support our observation of the dominance of the low-frequency phonon modes in contributing to the Na⁺ ion mobility in the lattice. High participation ratio values near the low-frequency phonon branches indicate that the movement of the Na⁺ ions is delocalized in those bands, which means that those vibrations are collective motions of many atoms, including Na⁺ ions and their neighbors (host sublattice), an easy way for ions to migrate, as the local energy landscape will change significantly when many atoms participate in the vibration. On the other hand, high-frequency phonon modes with low participation ratios suggest that the vibrations of Na⁺ ions are localized. Such localized vibrations generally correspond to small random displacements of Na⁺ ions (left panel in Fig. 4) and thus do not facilitate ion migration.

Fig. 4 (left panel) Phonon dispersions along high symmetry paths, (middle panel) frequency dependent phonon participation ratio, and (right panel) Na⁺ ion transport channel represented in red for the four representative structures with high diffusion coefficient materials: (a)–(c) Na₄TeS (OQMD ID 1554824), (d)–(f) NaVPdS₄ (OQMD ID 1391859), (g)–(i) Na₃YBr₆ (OQMD ID 1750683), and (j)–(l) Na₃ScBr₆ (ICSD ID 401335). The dispersions and participation ratio are colored by the magnitude of the mean squared displacements of Na⁺ ions contributed by each individual phonon mode.

Furthermore, the Na⁺ ion transport channels within the four representative structures are illustrated in the right panel of Fig. 4 (red isosurfaces), showing 3D diffusion in Fig. 4(c), (i), and (l), and 1D diffusion in Fig. 4(f). The red region is calculated based on the probability density of Na⁺ ion positions along the diffusion pathways obtained from AIMD trajectory data. These isosurfaces indicate the spatial regions with high Na⁺ ion occupancy, outlining the preferred migration pathways through the crystal lattice. The host structure framework composed of S and Te (yellow and blue, respectively) in Fig. 4(c), S, Pd, and V (yellow, blue, and orange, respectively) in Fig. 4(f), Y and Br (yellow and blue, respectively) in Fig. 4(i), and Sc and Br (yellow and blue, respectively) in Fig. 4(l) remains relatively static compared to the mobile Na⁺ ions and defines the bottlenecks and percolation geometry that govern the ion transport dynamics.⁶¹ The interconnectivity and continuity of the red isosurfaces reflect the dimensional characteristics of the Na⁺ ion conduction. In structures exhibiting 3D diffusion (Fig. 4(c), (i) and (l)), the isosurfaces form a well-connected and spatially extended network, with Na⁺ ion migrations along multiple crystallographic directions in those materials. In contrast, in 1D diffusion (Fig. 4(f)), the isosurfaces align along a single axis, which indicates highly anisotropic transport behavior. These transport channels highlight the influence of structural geometry on Na⁺ ion migration.

In the left panel of Fig. 5, we present the coupling factor between Na⁺ ions and the host sublattice for the four representative structures. The red region (coupling factor close to 1) indicates in-phase coupling, and the violet region indicates out-of-phase coupling. Most of these couplings occur at the low-frequency phonon bands up to ∼1.5 THz, which suggests vibrational coherence of Na⁺ ions with the surrounding lattice, facilitating framework-assisted ionic diffusion⁵² through collective lattice dynamics participation. The green and cyan regions (coupling factor close to 0) represent uncorrelated vibrations between Na⁺ ions and the host sublattice. In the right panel of Fig. 5, we show the partial VDOS of mobile Na⁺ ions and the host sublattice. All four structures exhibit large matched VDOS between Na⁺ ions and the sublattice in the broad frequency range. Such a feature is desirable for fast ionic transport. An interesting case is the material Na₄TeS (OQMD ID 1554824), where strong coupling between Na⁺ ions and the host sublattice can occur up to ∼2.5 THz and is also found in high frequency regions (above 3 THz). This can be understood from the partial VDOS shown in Fig. 5(b), where the VDOS between Na⁺ ions and the sublattice is matched over the entire frequency range. In particular, the center frequency of the partial VDOS of Na⁺ ions (∼3 THz) is about 2 times that of the host sublattice (∼1.4 THz). All 4 materials shown in Fig. 5 possess a high ratio of VDOS center frequency of Na⁺ ions to the acoustic cutoff frequency of the host sublattice, which is favorable for ionic transport as found earlier in Fig. 3.

Fig. 5 (left panel) Phonon dispersions colored by a Na⁺–sublattice coupling factor for each phonon mode. (right panel) Partial vibrational density of states (VDOS) with acoustic cutoff frequency of the structure and VDOS center of Na⁺ ions represented by black and magenta dashed lines, respectively. The host sublattice is denoted as “Non-Na⁺”. (a) and (b) Na₄TeS (OQMD ID 1554824), (c) and (d) NaVPdS₄ (OQMD ID 1391859), (e) and (f) Na₃YBr₆ (OQMD ID 1750683), and (g) and (h) Na₃ScBr₆ (ICSD ID 401335). Coupling factor values close to +1 (−1) indicate strongly correlated in-phase (out-of-phase) vibration between Na⁺ ions and host sublattices, while values close to 0 indicate uncorrelated vibration.

Before closing, we further calculated the elastic constants of the four representative materials by DFT, and their mechanical properties are presented in Table 2. The mechanical properties of the four superionic conductors, such as bulk modulus, shear modulus, Young's modulus, and hardness, are much lower than typical inorganic materials such as semiconductors and insulators^64–67 which again proves the soft lattice therein. We also observe that similar conclusions about the importance of lattice dynamics descriptors were reached for a larger class of solid-state electrolytes by López et al.,⁶⁸ who also highlighted a positive correlation between diffusion and vibrational descriptors, particularly those that explicitly incorporate anharmonic effects. This fact strongly vouches for the generality of the conclusions of the present study. While the present study focuses on validating the proposed lattice dynamics descriptors through high-throughput DFT- and AIMD-based benchmarking, a natural next step will be extending this framework toward screening of large-scale hypothetical structures to accelerate the discovery of new Na-superionic conductors.

Table 2 Mechanical properties of the four representative structures evaluated using DFT

Database	Material ID	Formula	Bulk modulus (GPa)	Shear modulus (GPa)	Young's modulus (GPa)	Poisson's ratio	Hardness (GPa)
OQMD	1554824	Na4TeS	22.93	13.92	34.74	0.25	3.37
OQMD	1391859	NaVPdS4	19.11	7.86	20.74	0.32	1.44
OQMD	1750683	Na3YBr6	14.32	5.89	15.55	0.32	1.18
ICSD	401335	Na3ScBr6	15.23	6.47	17.00	0.31	1.30

---

4 Conclusions

In summary, we establish a quantitative and strong correlation between the lattice dynamics and the diffusion coefficients based on universal machine learning potential equipped MD simulations and phonon calculations of 921 sodium structures, with some selected materials confirmed by direct DFT calculations and AIMD simulations. Large mean squared displacements of the material support faster Na⁺ ion transport, which is desirable for solid-state electrolytes in sodium all-solid-state battery applications. Additionally, we presented the significance of other phonon features, such as the acoustic cutoff frequency and the center phonon partial density of state of Na⁺ ions, in impacting the Na⁺ ion diffusion process. Lower acoustic cutoff frequency and lower center PDOS of Na⁺ ions just above the acoustic cutoff frequency correlate with high Na⁺ ion diffusion, which indicates that the lattice is softer and better matched vibrations between moving Na⁺ ions and the hosting sublattice in the low frequency range facilitate free motion of Na⁺ ions within the lattice, reinforcing the link between lattice dynamics and ionic mobility. We finally conducted phonon mode level analysis of the contribution to MSDs of Na⁺ ions and defined a coupling factor to quantitatively characterize correlated in-phase and out-of-phase vibrations between Na⁺ ions and host sublattices. Only extremely low frequency acoustic phonon modes and a limited number of low-lying low-energy optic phonon modes with strong Na⁺–sublattice coupling are primarily responsible for the large MSDs of Na⁺ ions, while the majority of phonon modes have little contribution. These insights from a lattice dynamics point of view highlight important features that can aid the design and accelerate the discovery of novel sodium superionic conductors with enhanced ionic conductivity, such as implementing the lattice dynamics features into machine learning-assisted material screening and inverse material design.

Author contributions

M. H. conveyed the idea and designed and supervised the study. O. A. performed MD simulations and processed data, and calculated lattice dynamics features. R. R. conducted significant DFT calculations for dynamic stability validation. M. A. developed the codes for running MD simulations with EquiformerV2 machine learning potential and the codes for phonon mode level analysis of MSD contributions. J. O. participated in postprocessing MD simulations and IFC data and Pearson correlation analysis. D. J. provided fruitful discussions and constructive suggestions on result analysis and presentation. O. A. prepared the draft of the manuscript. M. H., R. R., and D. J. revised the manuscript. All authors contributed to discussions and interpretation of results in the manuscript.

Conflicts of interest

The authors declare no competing financial or non-financial interests.

Data availability

All data reported herein are available upon reasonable request to the corresponding author.

Supplementary information is available. See DOI: https://doi.org/10.1039/d5mh01176k.

Acknowledgements

This work was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award # DE-SC0023377. Calculations were performed at Theia, an AI-focused HPC cluster at the University of South Carolina supported by NSF (award number 2320292), and at the Centro de Supercomputación de Galicia (CESGA) and at the Barcelona Supercomputing Center (BSC) within actions FI-2025-1-0010 and FI-2025-2-0015 of the Red Española de Supercomputación (RES).

References

1. M. Kim, K. Oka, S. Ahmed and D. A. Keen, Disordering Phenomena in Superionic Conductors Evidence for Superionic H₂O and Diffusive He–H₂O at High Temperature and High Pressure Disordering Phenomena in Superionic Conductors, J. Phys.: Condens. Matter, 2022, 34(39), 394001.
2. A. K. Sagotra, D. Chu and C. Cazorla, Influence of lattice dynamics on lithium-ion conductivity: A first-principles study, Phys. Rev. Mater., 2019, 3(3), 035405.
3. S. Wang, et al., Design principles for sodium superionic conductors, Nat. Commun., 2023, 14(1), 7615.
4. H. Kwak, et al., Emerging Halide Superionic Conductors for All-Solid-State Batteries: Design, Synthesis, and Practical Applications, ACS Energy Lett., 2022, 7, 1776–1805,  DOI:10.1021/acsenergylett.2c00438.
5. T. P. Bailey and C. Uher, Potential for superionic conductors in thermoelectric applications, Curr. Opin. Green Sustainable Chem., 2017, 4, 58–63,  DOI:10.1016/j.cogsc.2017.02.007.
6. R. M. Ormerod, Solid oxide fuel cells, Chem. Soc. Rev., 2003, 32, 17–28,  10.1039/b105764m.
7. A. J. E. Rettie, et al., A two-dimensional type I superionic conductor, Nat. Mater., 2021, 20, 1683–1688.
8. J. Kim, S. Kang and K. Min, Screening Platform for Promising Na Superionic Conductors for Na-Ion Solid-State Electrolytes, ACS Appl. Mater. Interfaces, 2023, 15, 41417–41425.
9. H. Gao, S. Xin, L. Xue and J. B. Goodenough, Stabilizing a High-Energy-Density Rechargeable Sodium Battery with a Solid Electrolyte, Chem, 2018, 4, 833–844.
10. PbSO, P., Koh, N., -, -e, Li Si, V. & LiX, S. PbO₂H₂OH₂SO₄/H₂OH + PbSO₄SO₄²⁻ + -e-2.0 V Lead-Acid Li_1−xCoO₂ LiPF₆ EC-DMC LiPF₆ EC-DMC Lithium Ion Future Batteries? Lithium Ion Based on Nanomaterials Ionic Liquids Lithium Metal Lithium-Air Lithium-Based Electrolytes, 1859.
11. X. He, Y. Zhu and Y. Mo, Origin of fast ion diffusion in super-ionic conductors, Nat. Commun., 2017, 8(1), 15893.
12. F. A. L. Laskowski, D. B. McHaffie and K. A. See, Identification of potential solid-state Li-ion conductors with semi-supervised learning, Energy Environ. Sci., 2023, 16, 1264–1276.
13. W. D. Richards, et al., Design and synthesis of the superionic conductor Na₁₀SnP₂S₁₂, Nat. Commun., 2016, 7(1), 11009.
14. Y. Wang, W. D. Richards, S. H. Bo, L. J. Miara and G. Ceder, Computational Prediction and Evaluation of Solid-State Sodium Superionic Conductors Na₇P₃X11 (X = O, S, Se), Chem. Mater., 2017, 29, 7475–7482.
15. Y. Xu, Y. Zong and K. Hippalgaonkar, Machine learning-assisted cross-domain prediction of ionic conductivity in sodium and lithium-based superionic conductors using facile descriptors, J. Phys. Commun., 2020, 4(5), 055015.
16. D. Park, et al., Computational screening of sodium solid electrolytes through unsupervised learning, npj Comput. Mater., 2024, 10(1), 226.
17. E. D. Cubuk, A. D. Sendek and E. J. Reed, Screening billions of candidates for solid lithium-ion conductors: A transfer learning approach for small data, J. Chem. Phys., 2019, 150(21), 214701.
18. F. A. L. Laskowski, D. B. McHaffie and K. A. See, Identification of potential solid-state Li-ion conductors with semi-supervised learning, Energy Environ. Sci., 2023, 16, 1264–1276.
19. A. D. Sendek, et al., Holistic computational structure screening of more than 12 000 candidates for solid lithium-ion conductor materials, Energy Environ. Sci., 2017, 10, 306–320.
20. Y. Liu, Q. Zhou and G. Cui, Machine Learning Boosting the Development of Advanced Lithium Batteries, Small Methods, 2021, 5(8), 2100442,  DOI:10.1002/smtd.202100442.
21. K. Fujimura, et al., Accelerated materials design of lithium superionic conductors based on first-principles calculations and machine learning algorithms, Adv. Energy Mater., 2013, 3, 980–985.
22. X. Wang, R. Xiao, H. Li and L. Chen, Oxysulfide LiAlSO: A Lithium Superionic Conductor from First Principles, Phys. Rev. Lett., 2017, 118(19), 195901.
23. Z. Yu, et al., A quaternary sodium superionic conductor – Na_10.8Sn_1.9PS1_1.8, Nano Energy, 2018, 47, 325–330.
24. Q. Zhou, et al., Sodium Superionic Conductors (NASICONs) as Cathode Materials for Sodium-Ion Batteries, Electrochem. Energy Rev., 2021, 4, 793–823,  DOI:10.1007/s41918-021-00120-8.
25. W. Wang, B. Jiang, L. Hu and S. Jiao, Nasicon material NaZr₂(PO₄)₃: A novel storage material for sodium-ion batteries, J. Mater. Chem. A, 2014, 2, 1341–1345.
26. A. M. A. Fami, N. A. Wahab, M. S. A. Rani, M. K. Yaakob and N. A. Mustaffa, First Principles Investigation of NASICON-Structured LiTi₂(PO₄)₃ and Mg_0.5Ti₂(PO₄)₃ Solid Electrolytes, Int. J. Electrochem. Sci., 2022, 17(1), 220115.
27. Y. Xu, Y. Zong and K. Hippalgaonkar, Machine learning-assisted cross-domain prediction of ionic conductivity in sodium and lithium-based superionic conductors using facile descriptors, J. Phys. Commun., 2020, 4(5), 055015.
28. S. Kang, M. Kim and K. Min Machine Learning-Aided Discovery of Superionic Solid-State Electrolyte for Li-Ion Batteries.
29. J. Lee, A. Seko, K. Shitara, K. Nakayama and I. Tanaka, Prediction model of band gap for inorganic compounds by combination of density functional theory calculations and machine learning techniques, Phys. Rev. B, 2016, 93(11), 115104.
30. S. Kirklin, et al., The Open Quantum Materials Database (OQMD): assessing the accuracy of DFT formation energies, npj Comput. Mater., 2015, 1, 15010.
31. A. Belsky, M. Hellenbrandt, V. L. Karen and P. Luksch, New developments in the Inorganic Crystal Structure Database (ICSD): accessibility in support of materials research and design, Acta Crystallogr., Sect. B: Struct. Sci., 2002, 58, 364–369.
32. L. Barroso-Luque, et al., Open Materials 2024 (OMat24) Inorganic Materials Dataset and Models, 2024.
33. J. Ojih, et al., High-throughput computational discovery of 3218 ultralow thermal conductivity and dynamically stable materials by dual machine learning models, J. Mater. Chem. A, 2023, 11, 24169–24183.
34. A. Rodriguez, et al., Unlocking phonon properties of a large and diverse set of cubic crystals by indirect bottom-up machine learning approach, Commun. Mater., 2023, 4, 61.
35. A. Rodriguez, et al., Million-scale data integrated deep neural network for phonon properties of heuslers spanning the periodic table, npj Comput. Mater., 2023, 9, 20.
36. M. Al-Fahdi, C. Lin, C. Shen, H. Zhang and M. Hu, Rapid prediction of phonon density of states by crystal attention graph neural network and high-throughput screening of candidate substrates for wide bandgap electronic cooling, Mater. Today Phys., 2025, 50, 101632.
37. M. Al-Fahdi, K. Yuan, Y. Yao, R. Rurali and M. Hu, High-throughput thermoelectric materials screening by deep convolutional neural network with fused orbital field matrix and composition descriptors, Appl. Phys. Rev., 2024, 11(2), 021402.
38. L. Chanussot, et al., Open Catalyst 2020 (OC20) Dataset and Community Challenges, ACS Catal., 2021, 11, 6059–6072.
39. H. Yang, et al., MatterSim: A Deep Learning Atomistic Model Across Elements, Temperatures and Pressures, 2024.
40. I. Batatia, D. P. Kovács, G. N. C. Simm, C. Ortner and G. Csányi MACE: Higher Order Equivariant Message Passing Neural Networks for Fast and Accurate Force Fields.
41. B. Deng, et al., CHGNet as a pretrained universal neural network potential for charge-informed atomistic modelling, Nat. Mach. Intell., 2023, 5, 1031–1041.
42. A. Togo and I. Tanaka, First principles phonon calculations in materials science, Scr. Mater., 2015, 108, 1–5.
43. S. Takamoto, et al., Towards universal neural network potential for material discovery applicable to arbitrary combination of 45 elements, Nat. Commun., 2022, 13, 2991.
44. R. Jacobs, et al., A practical guide to machine learning interatomic potentials – Status and future, Curr. Opin. Solid State Mater. Sci., 2025, 35, 101214.
45. F. Zhou, W. Nielson, Y. Xia and V. Ozoliņš, Compressive sensing lattice dynamics. I. General formalism, Phys. Rev. B, 2019, 100(18), 184308.
46. F. Zhou, W. Nielson, Y. Xia and V. Ozoliņš, Lattice anharmonicity and thermal conductivity from compressive sensing of first-principles calculations, Phys. Rev. Lett., 2014, 113(18), 185501.
47. F. Zhou, B. Sadigh, D. Åberg, Y. Xia and V. Ozoliņš, Compressive sensing lattice dynamics. II. Efficient phonon calculations and long-range interactions, Phys. Rev. B, 2019, 100(18), 184309.
48. C. Lin, J. Han, B. Xu and N. Marzari, First-principles phonon physics using the Pheasy code, 2025.
49. T. Krauskopf, et al., Comparing the Descriptors for Investigating the Influence of Lattice Dynamics on Ionic Transport Using the Superionic Conductor Na₃PS_4−xSe_x, J. Am. Chem. Soc., 2018, 140, 14464–14473.
50. S. Muy, et al., Tuning mobility and stability of lithium ion conductors based on lattice dynamics, Energy Environ. Sci., 2018, 11, 850–859.
51. C. Lin, X. Chen and X. Zou, Phonon–Grain-Boundary-Interaction-Mediated Thermal Transport in Two-Dimensional Polycrystalline MoS₂, ACS Appl. Mater. Interfaces, 2019, 11, 25547–25555.
52. M. K. Gupta, et al., Strongly Anharmonic Phonons and Their Role in Superionic Diffusion and Ultralow Thermal Conductivity of Cu₇PSe₆, Adv. Energy Mater., 2022, 12(23), 2200596.
53. T. Krauskopf, C. Pompe, M. A. Kraft and W. G. Zeier, Influence of Lattice Dynamics on Na⁺ Transport in the Solid Electrolyte Na₃PS_4−xSe_x, Chem. Mater., 2017, 29, 8859–8869.
54. M. K. Gupta, et al., Investigation of Low-Energy Lattice Dynamics and Their Role in Superionic Na Diffusion and Ultralow Thermal Conductivity of Na₃PSe₄ as a Solid-State Electrolyte, Chem. Mater., 2024, 36(23), 11377–11392.
55. T. Famprikis, et al., A New Superionic Plastic Polymorph of the Na⁺ Conductor Na₃PS₄, ACS Mater. Lett., 2019, 1, 641–646.
56. S.-H. Bo, Y. Wang, J. C. Kim, W. D. Richards and G. Ceder, Computational and Experimental Investigations of Na-Ion Conduction in Cubic Na₃PSe₄, Chem. Mater., 2016, 28, 252–258.
57. S. Xiong, et al., Na₃SbSe_4−xS_x as Sodium Superionic Conductors, Sci. Rep., 2018, 8, 9146.
58. N. A. S. Sampaio, et al., Applications of Correlation Analysis in Environmental Problems, Rev. Gest. Ambient. Sustentabilidade, 2024, 18, e04925.
59. D. Vivona, et al., Predicting Fast Oxygen Ion Conductors Via Descriptors for Migration Barrier Based on Electronic Structure and Lattice Dynamics. ECS Meeting Abstracts MA2023-01, 15–15 (2023).
60. S. Muy, et al., High-Throughput Screening of Solid-State Li-Ion Conductors Using Lattice-Dynamics Descriptors, iScience, 2019, 16, 270–282.
61. K. Gordiz, S. Muy, W. G. Zeier, Y. Shao-Horn and A. Henry, Enhancement of ion diffusion by targeted phonon excitation, Cell Rep. Phys. Sci., 2021, 2(5), 100431.
62. M. K. Gupta, et al., Fast Na diffusion and anharmonic phonon dynamics in superionic Na₃PS₄, Energy Environ. Sci., 2021, 14, 6554–6563.
63. M. A. Kraft, et al., Influence of Lattice Polarizability on the Ionic Conductivity in the Lithium Superionic Argyrodites Li₆PS₅X (X = Cl, Br, I), J. Am. Chem. Soc., 2017, 139, 10909–10918.
64. S. Chibani and F.-X. Coudert, Systematic exploration of the mechanical properties of 13 621 inorganic compounds, Chem. Sci., 2019, 10, 8589–8599.
65. M. de Jong, et al., Charting the complete elastic properties of inorganic crystalline compounds, Sci. Data, 2015, 2, 150009.
66. E. Kiely, R. Zwane, R. Fox, A. M. Reilly and S. Guerin, Density functional theory predictions of the mechanical properties of crystalline materials, CrystEngComm, 2021, 23, 5697–5710.
67. K. Li, P. Yang, L. Niu and D. Xue, Hardness of Inorganic Functional Materials, Rev. Adv. Sci. Eng., 2012, 1, 265–279.
68. C. López, A. Emperador, E. Saucedo, R. Rurali and C. Cazorla, Universal ion-transport descriptors and classes of inorganic solid-state electrolytes, Mater. Horiz., 2023, 10, 1757–1768.
69. Z. Zhang, et al., Na₁₁Sn₂PS₁₂: A new solid state sodium superionic conductor, Energy Environ. Sci., 2018, 11, 87–93.
70. A. Banerjee, et al., Na₃SbS₄: A Solution Processable Sodium Superionic Conductor for All-Solid-State Sodium-Ion Batteries, Angew. Chem., Int. Ed., 2016, 55, 9634–9638.
71. S. Takeuchi, K. Suzuki, M. Hirayama and R. Kanno, Sodium superionic conduction in tetragonal Na₃PS₄, J. Solid State Chem., 2018, 265, 353–358.

---