Chemical similarity description of liquid kinetics across composition space in magnesium alloys

Ab initio molecular dynamics simulations

Binary Mg–X liquid alloys were systematically investigated by ab initio molecular dynamics (AIMD) simulations, where X denotes 31 common alloying elements in Mg alloys, including Li, Be, Al, Si, Ca, Sc, Mn, Fe, Ni, Cu, Zn, Ga, Sr, Y, Zr, Ag, Sn, Sb, Bi, La, Ce, Pr, Nd, Sm, Gd, Tb, Dy, Ho, Er, Yb, and Th. A statistical analysis of Mg alloy compositions reported in the literature [Supplementary Figure 1] indicates that the solute concentration is typically below 10 at.%. Accordingly, two solute concentrations, 2 and 5 at.%, were considered in this study. The simulated temperatures were selected as 1,500, 1,200, and 1,000 K, all above the melting point of pure Mg (923 K), to ensure a liquid state.

The AIMD simulations were performed using the Vienna ab initio simulation package (VASP)^[18] within the generalized gradient approximation (GGA) and the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional^[19]. Core electrons were replaced by the projector augmented wave (PAW) pseudopotentials^[20], treating Mg (3s) as the valence state. Spin-polarized calculations were used for transition and rare-earth metal elements. The specific PAW potentials employed for each solute element are listed in Table 1. Their selection followed the recommendations of the official VASP pseudopotential documentation while considering the computational cost of the extensive AIMD simulations. The energy cutoff was set to 300 eV, and a single Γ point was adopted to sample the Brillouin zone. AIMD simulations were carried out in the canonical ensemble (NVT), i.e., at a constant number of particles (N), constant volume (V), and a constant temperature (T), using a cubic supercell containing 200 atoms, which was equilibrated for 20 ps at the target temperature with a timestep of 2 fs. This NVT trajectory was used for analysis, and no additional production simulation was performed. In total, 186 AIMD simulations (31 systems × 2 concentrations × 3 temperatures) were conducted, covering a broad composition space, with a total simulation time of 3.72 ns. The data analysis was performed using the PyMatterSim code^[21]. All figures were prepared using Matplotlib^[22] and Origin Learning Edition.

Calculation of liquid kinetics

To quantify the kinetic properties of the liquid alloys, the self-diffusion coefficient D was derived from the calculated mean squared displacement (MSD) according to the Einstein relation^[23],

where r is the atomic position vector, t is time, N represents the number of atoms, and the angular brackets denote ensemble averaging. Note that to improve statistical reliability, MSDs were computed using a moving average with a 2 ps window width over 20 ps trajectories for each system. The diffusion coefficient D can be estimated from the slope of the MSD in the long-time limit:

The statistical uncertainty was quantified as the standard deviation of the D obtained from the individual moving windows. The diffusion coefficient uncertainty was quantified as the standard deviation over the moving windows. We calculated the average self-diffusion coefficient of the liquid alloy, obtained from the total MSD of all atoms and thus characterizing the overall atomic mobility of the system.

Symbolic regression

Symbolic regression is a data-driven machine learning method that aims to construct explicit analytical expressions that best describe the relationship between input variables and target properties, without assuming a predefined functional form^[24]. In this work, we implemented PySR code^[25] to construct a symbolic regression model for predicting the calculated diffusion coefficients by elemental properties, including atomic mass, radius, period number, group number, valence number, electronegativity, and temperature. The solute concentration was incorporated as a weighting factor applied to the elemental descriptors, thereby enabling the model to account for composition-dependent effects. Candidate expressions were iteratively constructed and optimized using a set of mathematical operators - unary operators including “square”, “cube”, and “exp”, and binary operators including “+”, “-”, “×”, “/” - to achieve a balance between predictive accuracy and model simplicity.

High-throughput computational workflow

Determining the zero-pressure equilibrium volume across the broad composition space investigated is nontrivial, given that conventional approaches rely on trial-and-error adjustments of the cell size. To enable high-throughput calculations, we developed an automated workflow for equilibrium volume determination. As illustrated in Figure 1, the workflow begins with generating an initial configuration for a given alloy composition and temperature using an in-house script. Atomic positions were randomly distributed within a cubic simulation cell, while enforcing a minimum interatomic distance of 2.9 Å (approximately twice the atomic radius of Mg) to avoid unphysical overlaps. Initial velocities were sampled from the Maxwell-Boltzmann distribution corresponding to the target temperature. The system was first subjected to an isothermal–isobaric (NPT) simulation, i.e., at a constant number of particles (N), constant pressure (P), and a constant temperature (T), to obtain a physically reasonable initial cell volume, a suitable starting point for the subsequent zero-pressure optimization. Subsequently, NVT simulations were performed to iteratively refine the cell volume. The magnitude of each box-length adjustment is scaled with the number of atoms (N_atom) and empirically determined according to the pressure deviation. Specifically, pressure deviations of 50, 20, 7, and 3 correspond to adjustment coefficients of 8E-3, 4E-3, 1E-3, and 4E-4, respectively. If convergence (average pressure <P> within 2 kbar) was not achieved within three iterations, a quadratic pressure–volume relation (P = aV² + bV + c) was constructed by fitting the pressure–volume data from the NVT runs, and the equilibrium volume was determined at <P> = 0^[26]. The convergence of the AIMD simulations is illustrated in Supplementary Figure 2, which presents the time evolution of pressure and internal energy for the Mg–5 at.% Al liquid alloy at 1,500 K.

Figure 1. High-throughput computational workflow for determining zero-pressure equilibrium volume in AIMD simulations of Mg liquid alloys at given composition C and temperature T. AIMD: Ab initio molecular dynamics; NPT: isothermal–isobaric ensemble; NVT: canonical ensemble.

Liquid kinetics

Using the AIMD trajectories generated from the above computational workflow, we can systematically evaluate the diffusion coefficients D of binary Mg alloys. D serves as a fundamental kinetic parameter of liquid alloys, reflecting atomic mobility and governing mass transport phenomena such as diffusion-controlled mixing and segregation. We calculated D values from the MSD curves [Supplementary Figure 3], which exhibit well-defined linear regimes at long times, indicative of diffusive motion. To assess the reliability of the adopted computational setup, we first calculated the diffusion coefficient of pure liquid Mg. The calculated value is compared with available experimental data and previously reported AIMD and classical molecular dynamics (MD) results shown in Supplementary Figure 4. Our calculated diffusivity falls within the range of reported values, despite the differences in simulation and experimental temperatures, and is consistent with previous computational results.

A primary concern in calculating D is the potential for finite-size artifacts, where the periodic boundary conditions and limited cell dimensions may constrain long-range hydrodynamic correlations. To address this, we perform large-scale molecular dynamics simulations powered by deep learning interatomic potentials [deep potential (DP)] to evaluate finite-size effects. We have demonstrated that the DP model faithfully reproduces the liquid kinetics of the Mg–Al system in our previous work^[8], showing excellent agreement with both AIMD benchmarks and experimental measurements. See more details about DP development in the Supplementary Materials. The deep potential molecular dynamics (DPMD) simulations were initialized by replicating the AIMD-equilibrated supercells into 2 × 2 × 2 (1,600 atoms), 3 × 3 × 3 (5,400 atoms), and 4 × 4 × 4 (12,800 atoms) configurations. As illustrated in Figure 2, the scaled D/D_N=12,800 approaches unity across the investigated temperatures. For the AIMD simulation cell containing 200 atoms, the calculated diffusion coefficients deviate by less than 10% from those of the largest supercell. Specifically, the diffusion coefficients at 1,500, 1,200, and 1,000 K are 1.611, 0.997, and 0.624 Ų/ps for the 200-atom cell and 1.576, 0.935, and 0.639 Ų/ps for the 12,800-atom cell, with relative deviations of 2.2%, 6.2%, and 2.3%, respectively, indicating that finite-size effects are small for the Mg–5 at.% Al system under the present simulation conditions. Although chemistry-dependent finite-size effects cannot be completely excluded for other Mg–X systems, this result supports the use of 200-atom cells for the AIMD simulations in this study. It should also be acknowledged that the diffusion behavior presented herein focuses on the total diffusion coefficients.

Figure 2. Finite-size effect in determining the diffusion coefficient of Mg–5 at.% Al liquid alloys at temperatures 1,500, 1,200, and 1,000 K.

The temperature dependence of the diffusion coefficient follows the Arrhenius law, D = D₀exp(-Q/RT), where D₀ and Q are the pre-exponential factor and activation energy, respectively. Note that D and D₀ are reported in Ų/ps throughout this work, consistent with the characteristic length and time scales of the MD simulations. For comparison with experimental and literature data, these quantities can be converted to the conventional unit of cm²/s using 1 Ų/ps = 10^-4 cm²/s. Arrhenius fits of D for 31 binary Mg-based liquid alloys are illustrated in Figure 3, with the corresponding D₀ and Q summarized in Table 1. With increasing solute concentration from 2 to 5 at.%, the D values decrease for most systems [Supplementary Figure 5]. Despite dilute solute concentrations, binary Mg-based liquid alloys exhibit distinct diffusivity variations across different solutes, indicating a strong dependence of liquid kinetics on alloying chemistry. A natural question is whether such variations can be rationalized by elemental descriptors. From a classical kinetic-theory perspective, the diffusion coefficient of metallic liquids can be approximated as inversely proportional to the square root of atomic mass $$ D\propto \frac{1}{\sqrt{m}} $$, where m is the atomic mass^[14]. This trend is partially consistent with the relatively high diffusivity observed in binary liquids containing light solutes, such as Mg–Li and Mg–Be. However, this simple scaling fails to capture the overall behavior. Notably, binary liquids containing heavier solutes such as Sn, Sb, and Bi exhibit unexpectedly high diffusivity, comparable to or even exceeding that of the light-solute systems. Such deviations indicate that atomic mass alone is insufficient to describe diffusion in Mg-based liquid alloys. As shown in Supplementary Figure 6, the diffusion coefficients exhibit no obvious correlation with either atomic radius or atomic weight for binary Mg liquid alloys. These results highlight the need for a chemistry-informed descriptor capable of capturing the dependence of diffusivity on alloying chemistry.

Figure 3. Arrhenius fits of the diffusion coefficients for binary Mg liquid alloys at solute concentrations (A) 2 at.% and (B) 5 at.%.

Chemical similarity description of liquid kinetics

To elucidate the dependence of liquid kinetics on chemistry, we first examine the periodic-table mapping of the calculated diffusion coefficients for binary Mg liquid alloys, as illustrated in Figure 4. A clear organization emerges, wherein elements located within specific regions of the periodic table exhibit similar liquid kinetics. In particular, binary liquid alloys with s-block solutes (alkali and alkaline-earth metals) and early p-block transition metal solutes exhibit relatively high D (marked by red in Figure 4). In contrast, the binary systems with d-block or f-block rare-earth elements show markedly lower D (marked by blue in Figure 4). We acknowledge that this classification reflects an overall trend rather than a universal rule, as deviations occur for some elements. The observed periodic organization of liquid kinetics may originate from the underlying chemical similarities among solute elements across the periodic table. Similarities in electronic configuration, atomic size, and electronegativity among neighboring elements in the periodic table promote analogous interactions with the Mg solvent, resulting in comparable effects on liquid-state atomic transport.

Figure 4. Periodic-table mapping of the diffusion coefficients for binary Mg liquid alloys at solute concentrations of 2 at.% (top row) and 5 at.% (bottom row) and temperatures of 1,500, 1,200, and 1,000 K. The color scale ranges from blue (low diffusivity) to red (high diffusivity).

While this representation provides an intuitive qualitative overview, quantitative description and predictive modeling require a continuous descriptor that captures the underlying chemical similarities among solute elements. To this end, we introduce the Pettifor chemical scale^[27-29], a physically motivated yet phenomenological ordering in which chemically related elements occupy neighboring positions. It was originally proposed to rationalize the structural stability of crystalline compounds^[27,28], and is built on the empirical observation that chemical similarity governs preferences for crystalline prototypes. The Pettifor chemical scale values of the solute elements are listed in Table 1. As shown in Figure 5, projecting the two-dimensional periodic trends onto this one-dimensional scale reveals a remarkably simple behavior: the diffusion coefficient depends approximately quadratically on the Pettifor chemical scale across all investigated temperatures and solute concentrations. The high coefficient of determination (R²) attests to the robustness of this relationship, indicating that the complex chemical dependence of liquid kinetics can be largely captured by this reduced descriptor. This finding extends the Pettifor chemical scale from static crystalline structural stability to the dynamic behavior of liquid alloys, suggesting a unified chemical origin underlying both thermodynamic and dynamic properties.

Figure 5. Quadratic scaling of the diffusion coefficient with the Pettifor chemical scale of the solute element for corresponding binary Mg liquid alloys, across the investigated concentrations and temperatures. Solid lines indicate second-order polynomial fits, demonstrating a non-monotonic dependence. R²: Coefficient of determination.

The observed quadratic dependence on the Pettifor chemical scale is an empirical relationship that organizes the liquid diffusivities of the Mg–X systems examined in this work. While this smooth variation may ultimately reflect systematic changes in solute–solvent interactions along the chemical coordinate, the present analysis does not directly establish the underlying physical mechanisms. The chemical-similarity model is therefore best regarded as a phenomenological description that captures the dominant chemical tendency.

We first revisit the frameworks for investigating the dependence of liquid kinetics on alloying chemistry as shown in Figure 6A. Building on the empirical quadratic relationship with the Pettifor chemical scale, we construct a phenomenological model that simultaneously accounts for both compositional and temperature effects on liquid diffusivity, which we refer to as the chemical similarity description model

Figure 6. (A) Frameworks for investigating the dependence of liquid kinetics on alloying chemistry; The prediction models of diffusion coefficients across composition space by (B) our proposed chemical similarity description model, (C) data-driven symbolic regression machine learning model, and (D) physics-based excess entropy scaling. R²: Coefficient of determination; MAE: mean absolute error.

where P denotes the Pettifor chemical scale and C is the solute concentration. Note that the solute concentration C is expressed in at.%, i.e., C = 2 and 5 for 2 at.% and 5 at.% systems, respectively. Motivated by previous observations that the diffusivity of simple liquids approximately follows a power-law temperature dependence $$ D\propto T^n $$ (with n typically between 1.7 and 2.3)^[14], we adopt the form $$ (\frac{T}{T_0})^{\alpha} $$. The fitting parameters are a, b, c, α, and T₀ with values of 1.52 × 10^-4, 1.39 × 10^-2, 3.34 × 10^-2, 1.92, and 1,072, respectively. The model achieves an R² of 0.983 and a mean absolute error (MAE) of 0.044 on the dataset used for parameter fitting [Figure 6B]. In addition, our proposed description model also preserves clear physical interpretability. Both P and C are dimensionless quantities, and the temperature term $$ (\frac{T}{T_0})^{\alpha} $$ is also dimensionless by construction. The fitted reference temperature T₀ = 1,072 K is primarily used to normalize the power-law temperature dependence. Its proximity to the characteristic liquidus range of the studied Mg-based alloys may suggest a possible connection with the thermodynamic behavior of the liquid phase. However, since T₀ is obtained from a multivariable phenomenological fit, this interpretation cannot be established quantitatively here. We therefore treat T₀ as a fitted reference temperature.

To assess the predictive capability of the chemical similarity description model, we performed a hold-out validation test. Specifically, the data associated with Al, Zn, Ag, and Sm were excluded from the fitting process, and the model parameters were determined from the remaining Mg–X systems. The resulting model was then used to predict the diffusivities of these four held-out solutes across the investigated temperatures and concentrations. As shown in Supplementary Figure 7, the predicted diffusivities agree closely with the corresponding AIMD results (R² = 0.984 and MAE = 0.045), demonstrating that the chemical similarity description retains predictive capability for Mg–X systems not included in model parameterization. We note that incorporating a broader range of solute elements and concentrations can improve the generality of the fitted parameters.

To provide a data-driven benchmark, machine learning model^[30] can establish a direct relationship between alloying chemistry and liquid kinetics (see lower panel in Figure 6A). In this work, we implement a symbolic regression model^[31,32]. The candidate elemental descriptors include atomic mass, radius, valence, period number, group number, and temperature. Each elemental property X is represented by its concentration-weighted effective value according to the alloy composition: X_eff = (1 - C_frac)X_Mg + C_fracX_solute, where C_frac is the solute atomic fraction, taken as 0.02 and 0.05. X_Mg and X_solute are the properties of pure Mg and the corresponding solute element, respectively. The resulting expression obtained by symbolic regression is:

where r (in Å), m (in amu), p (dimensionless) represent elemental radius, mass, and period number, respectively. As illustrated in Figure 6C, this approach also achieves high fitting accuracy (R² = 0.978, MAE = 0.050) across the investigated temperatures and compositions. Nevertheless, the obtained expression remains largely empirical, offering limited physical insight and lacking strict dimensional consistency.

Furthermore, we review the conventional frameworks used to link alloying chemistry with liquid kinetics. As illustrated in the upper panel of Figure 6A, most previous studies adopt a structure-mediated pathway: alloying chemistry modifies the liquid structure, which in turn governs atomic transport^[33,34]. A well-established quantitative bridge in this framework is the excess-entropy scaling relation^[11]. In this physically grounded description, liquid kinetics is controlled by the short-range repulsive interactions, which can be approximated by binary hard-sphere collisions. The frequency of structural relaxation in the liquid dictates cage diffusion and is proportional to the number of accessible configurations in the system, leading to

where σ is the hard sphere diameter, Γ is the collision frequency, S is the excess entropy per atom. C and β are fitting parameters with values of 0.039 and 1.080, respectively, obtained from the present liquid kinetics dataset. Based on the Enskog theory of atomic transport, Γ = 4σ²g(σ)ρ$$ \sqrt{\frac{\pi k_B T}{m}} $$, where g(σ) is the pair correlation function g(r) evaluated at σ, which may be interpreted at the position of the first peak in g(r), m is atomic mass, ρ is number density, and k_B is Boltzmann constant. The excess entropy can be restricted to the two-body approximation: S₂ = -2πρ∫₀^∞[g(r)lng(r) - g(r) + 1]r²dr. The dependence of the pair excess entropy S₂ on temperature and composition is illustrated in Supplementary Figure 8. As shown in Figure 6D, the excess entropy scaling captures the overall trend of diffusivity across different alloy systems. However, the fitting accuracy remains moderate (R² = 0.898, MAE = 0.101), with noticeable deviations for the 5 at.% systems. These discrepancies likely arise from the increased complexity of local chemical environments at higher solute concentrations, indicating the limitation of the quasi-hard-sphere approximation^[12]. Moreover, the reliance on structural sampling entails significant atomistic computational cost, which limits its applicability for predictive modeling. The correlations between liquid kinetics and other commonly used structural descriptors (including Voronoi volume^[35,36], bond length^[37], and bond orientation order^[38,39]) are shown in Supplementary Figures 9-11 for the studied binary liquid alloys.

A quantitative comparison of the three approaches is summarized in Table 2. Compared with excess-entropy scaling, which is physically rigorous and dimensionally consistent but offers only moderate predictive accuracy, the present chemical-similarity model achieves substantially higher accuracy without requiring any explicit structural input. In contrast to symbolic regression, which offers high predictive accuracy at the expense of physical interpretability, our model preserves clear physical meaning. Consequently, the chemical similarity description provides an optimal balance between physical insight and predictive power, rendering it highly attractive for efficient modeling of complex liquid alloy systems.