Iron oxidation is a complex process involving critical atomistic events, such as atomic adsorption, diffusion, and surface reconstruction, understanding of which is significant for both surface science and coating technology. Atomistic simulation serves as an useful tool to investigate the processes, where description of interatomic interactions is required. However, selecting appropriate force field or interatomic potential is not only difficult, but also essential for getting accurate result. In this work, we present a detailed benchmark of reactive force fields (ReaxFFs) and universal machine learning interatomic potentials (uMLIPs) against density functional theory (DFT) calculations of oxygen adsorption on various α-iron surfaces, which is the first yet crucial step towards oxidation. The comparisons show the coverage-dependent performance and improvable accuracy of both ReaxFFs and uMLIPs at reproducing DFT results, with ReaxFFs outperforming uMLIPs. Subsequently, iron oxidation is simulated using ReaxFF and uMLIP. The results reveal the strong capability of ReaxFF and poor stability of uMLIP for describing reactive process, i.e., the formation of iron oxide. This may be attributed to the suitable functional form of ReaxFF for the description of bond changes. The insights presented here not only provide an example of benchmarking force field or interatomic potential for system of interest, but also highlight the applicability of ReaxFF and scopes of improvement of uMLIP.

Universal machine-learning interatomic potentials (uMLIPs) are emerging as foundation models for atomistic simulation, offering near-ab initio accuracy at far lower cost. Their safe, broad deployment is limited by the absence of reliable, general uncertainty estimates. We present a unified, scalable uncertainty metric, U, built from a heterogeneous ensemble that reuses existing pretrained MLIPs. Across diverse chemistries and structures, U strongly tracks true prediction errors and robustly ranks configuration-level risk. Using U, we perform uncertainty-aware distillation to train system-specific potentials with far fewer labels: for tungsten, we match full density-functional-theory (DFT) training using 4% of the DFT data; for MoNbTaW, a dataset distilled by U supports high-accuracy potential training. By filtering numerical label noise, the distilled models can in some cases exceed the accuracy of the MLIPs trained on DFT data. This framework provides a practical reliability monitor and guides data selection and fine-tuning, enabling cost-efficient, accurate, and safer deployment of foundation models.

Enhancing the hydrogen embrittlement (HE) resistance of alloys caters to the urgent needs of engineering safety and long-distance hydrogen transportation. Highly dense precipitates in the alloys act as H traps, however, some of them cannot strongly trap H thus failing to prevent its accumulation at the critical regions. Experimentally, it is challenging to expeditiously identify and generate phases causing strengthening and acting as strong H traps. Here, we demonstrate a computation-based design strategy to generate precipitates strongly trapping H. Based on the quantum machine learning Al-Sc-Cu potential, the optimal processing parameters of strong H trapping phases are determined, even though they are metastable in nature. Elemental mapping in electron microscope and atom probe tomography confirms the presence of Cu in Al₃Sc and its strong interaction with H. Hence, we envisage the proposed strategy will accelerate the design of HE-resistant microstructures of various technologically relevant materials via identification of desirable phases.

Cold-drawn pearlite wire is widely used in industry due to its exceptional high strength. Understanding the deformation mechanisms during the cold-drawing process of pearlite, particularly the deformation and decomposition of cementite, is of great significance. In this study, a bicrystal model tailored to lamellar structures is developed to calculate the elastic properties and stress concentration of pearlite. By analyzing slip activation in both ferrite and cementite, along with the yield strength, we reveal the significant influence of loading direction on pearlite deformability. Notably, the yield strength varies from 9.5 GPa to 17.0 GPa. Under specific loading conditions, plastic deformation is observed to initiate in cementite, challenging the conventional assumption that slip bands always originate in ferrite. Furthermore, factors that influence the plastic deformation of pearlite are discussed. A successive strengthening mechanism is proposed to explain the excellent deformability and high strength of pearlite after extensive deformation. This work introduces a novel method for directional loading of lamellar structures. The surprising finding that plastic deformation, without fracture, can initiate in cementite, might offer directions for developing other structural materials with extreme tensile strength and deformability.

Extended defects such as dislocation networks and general grain boundaries are ubiquitous in metals, and accurate modeling these extensive defects is crucial to elucidate their deformation mechanisms. However, existing machine learning interatomic potentials (MLIPs) often fall short in adequately describing these defects, as their large characteristic scales exceed the computational limits of first-principles calculations. To address this challenge, we present a computational framework combining a defect genome constructed via empirical interatomic potential-guided sampling, with an automated reconstruction technique that enables accurate first-principles modeling of general defects by converting atomic clusters into periodic configurations. The effectiveness of this approach was validated through simulations of nanoindentation, tensile deformation, and fracture in BCC tungsten. This framework enhances the modeling accuracy of extended defects in crystalline materials and provides a robust foundation for advancing MLIP development by leveraging defect genomes strategically.

This extensive review highlights the central role of classical molecular simulation in advancing hydrogen (H2) technologies. As the transition to a sustainable energy landscape is urgently needed, the optimization of H2 processes, spanning production, purification, transportation, storage, safety, and utilization is essential. To this end, accurate prediction of thermodynamic, transport, structural, and interfacial properties is important for overcoming engineering challenges across the entire H2 value chain. Experimental measurements, despite being the traditional way of obtaining these properties, can be limited by the distinctive nature of H2, harsh operating conditions, safety constraints, and extensive parameter spaces. Free from such limitations, classical molecular simulations, in the general frameworks of Monte Carlo and Molecular Dynamics, provide an optimal balance between computational efficiency and accuracy, bridging the gap between quantum mechanical calculations and macro-scale modeling. This review also systematically covers molecular simulation methods and force fields for computing key properties of H2 systems, such as phase and adsorption equilibria and transport coefficients. Beyond property prediction, we explore how molecular simulation reveals fundamental mechanisms governing hydrate formation and dissociation, membrane permeations, and H2 embrittlement. When possible, data from multiple sources are compared and critically assessed, while effort is put on evaluating the force fields used and methodological approaches followed in the literature. Finally, this review aims at identifying research gaps and future opportunities, emphasizing emerging approaches, such as molecular simulation in the era of artificial intelligence.

Understanding atomic hydrogen (H) diffusion in multi-principal element alloys (MPEAs) is crucial for enhancing hydrogen transport and storage technologies. However, the vast compositional space and complex chemical environments of MPEAs pose significant challenges. We develop highly accurate machine learning force field and neural network-driven kinetic Monte Carlo simulations to investigate H diffusion in body-centered cubic (BCC) MoNbTaW MPEAs. H diffusion exhibits super-Arrhenius behavior in MPEAs, dominated by the low percentile of the H solution energy spectrum. Robust analytical models are derived via machine learning symbolic regression to predict H diffusivity across general BCC MPEAs. Additionally, it is revealed that chemical short-range order (SRO) generally does not impact H diffusion in MoNbTaW MPEAs, except it enhances diffusion when H-favoring elements are present in low concentrations. These insights not only deepen our understanding of H diffusion dynamics in MPEAs but also guide the strategic development of advanced MPEAs for hydrogen-related applications by manipulating element type, composition, and SRO.

In this study, we explore the mechanisms underlying the exceptional intrinsic strength of face-centered cubic (FCC) Multi-Principal Element Alloys (MPEAs) using a multifaceted approach. Our methods integrate atomistic simulations, informed by both embedded-atom model and neural network potentials, with first-principles calculations, stochastic Peierls-Nabarro (PN) modeling, and symbolic machine learning. We identify a consistent, robust linear correlation between the strength of MPEAs and the standard deviation of the maximum stacking-fault restoring force (τ_max,sd) across various potentials. This finding is substantiated by comparing the experimental strengths of Cantor alloys’ subsystems and Ni_62.5V_37.5 against τ_max,sd values from high-throughput first-principle calculations. Our theoretical insights are derived from integrating the stochastic Peierls-Nabarro model with a shearable precipitation hardening framework, demonstrating that lattice distortion alone does not directly enhance intrinsic strength. Instead, τ_max,sd emerges as a critical determinant, capable of boosting the strength of MPEAs by up to tenfold. Our analysis reveals the critical role of the exponential form of the PN model in achieving substantial strength improvement by transforming the Gaussian-like distribution of τ_max into an exponential-like distribution of local Peierls stress. Additionally, using an advanced symbolic machine learning technique, the sure independence screening and sparsifying operator (SISSO) method, we derive interpretable relationships between MPEA strength, elastic properties, and τ_max statistics, offering new insights into the design and optimization of advanced MPEAs. These findings highlight that the nonlinear physics and atomic fluctuations characterizing MPEAs not only underpin their unconventional intrinsic strength but also contribute to other complex properties such as sluggish diffusion and cocktail effect.

Multi-principal element alloys (MPEAs) are renowned for their enhanced mechanical strength relative to their constituent metals, as evidenced by various experimental techniques such as tension/compression tests and instrumental indentation. Nevertheless, atomistic simulations sometimes produce conflicting results, casting doubt on the consistently superior mechanical properties of MPEAs. In this study, machine-learning interatomic potentials (MLIPs) with first-principles accuracy were developed for body-centered cubic refractory MoNbTaW MPEAs, enabling systematic atomistic simulations under various deformation scenarios. The new MLIPs are supported by a comprehensive dataset encompassing extensive defects, and the established embedded-atom model (EAM) potential was benchmarked against both this dataset and the new MLIP. Simulations covering diverse compositions confirm that both MLIPs and EAM accurately capture the critical strengthening mechanisms in MoNbTaW MPEAs. It is revealed that MPEAs generally exhibit superior mechanical strength compared to their constituent metals in macroscale specimens, primarily due to solid solution strengthening during dislocation motion. However, at the nanoscale—where plasticity is predominantly governed by dislocation nucleation and grain boundary deformation—the constituent metals may outperform MPEAs. A critical length scale is identified above which MPEAs demonstrate enhanced mechanical strength relative to their constituent elements; below this scale, the advantage diminishes, underscoring a significant size-dependent effect that is crucial for optimizing MPEA applications, particularly at the nanoscale.

Recent advances in machine learning, combined with the generation of extensive density functional theory (DFT) datasets, have enabled the development of universal machine learning interatomic potentials (uMLIPs). These models offer broad applicability across the periodic table, achieving first-principles accuracy at a fraction of the computational cost of traditional DFT calculations. In this study, we demonstrate that state-of-the-art pretrained uMLIPs can effectively replace DFT for accurately modeling complex defects in a wide range of metals and alloys. Our investigation spans diverse scenarios, including grain boundaries and general defects in pure metals, defects in high-entropy alloys, hydrogen-alloy interactions, and solute-defect interactions. Remarkably, the latest EquiformerV2 models achieve DFT-level accuracy on comprehensive defect datasets, with root mean square errors below 5 meV atom−1 for energies and 100 meV Å−1 for forces, outperforming specialized machine learning potentials such as moment tensor potential and atomic cluster expansion. We also present a systematic analysis of accuracy versus computational cost and explore uncertainty quantification for uMLIPs. A detailed case study of tungsten (W) demonstrates that data on pure W alone is insufficient for modeling complex defects in uMLIPs, underscoring the critical importance of advanced machine learning architectures and diverse datasets, which include over 100 million structures spanning all elements. These findings establish uMLIPs as a robust alternative to DFT and a transformative tool for accelerating the discovery and design of high-performance materials.

Prone H reduction is considered an important factor in the poor corrosion resistance of Mg and its alloys, while the reduced H simultaneously impacts their mechanical properties whose mechanism is still unclear. It can be experimentally found that the elongation of Mg charged with atomic H is 2.76 % greater than that in air. To reveal the underlying physics, multi-scale modeling combining first-principle calculation, molecular dynamic/static (MD/MS) simulation, and crystal plasticity finite element method (CPFEM) is first employed to elaborate the influence of H on Mg at different length scales. The first-principle results show that the Prism-I {101¯0} exhibits the most corrosive nature with an effective H adsorption density that reaches 18 nm⁻² and its diffusion barrier is only 0.156 eV H⁻¹. Conversely, the Basal {0001} has the best surficial H resistance. After H infiltration into the Mg matrix, the generalized stacking fault energies of most twining planes decrease by 2.26 % ∼18.49 %. Especially for the Basal {0001}, the H not only lowers its stacking fault energy to -7.13 J m⁻², but also impedes its cleavage cracking along [101¯0] according to the MD/MS simulation. The presence of H within the grains induces early initiation of stacking fault and elevates the critical stress at the crack tips. The CPFEM modeling reveals that the difference in twining growth is concentrated within 4 % strain. The H addition promotes the twining of Mg, however, following 4 % strain, the relative activity of planes in the Mg/Mg-H models is consistent.

In multiscale modeling methods (MMM), the integration of atomistic to continuum coupling is a common practice, where two regions have their own distinct length scales. The equilibrium configurations of such multiscale systems under given conditions are typically obtained through energy minimization algorithms (EMA). However, traditional EMAs, such as the conjugate gradient (CG) and limited-memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) algorithms, are unable to discern the diverse scales inherent in such systems. In this work, it is found that the convergence rate of energy minimization in multiscale simulations is significantly slower than that in full atomistic simulations, regardless of using CG, LBFGS algorithms or the latest fast inertial relaxation engine (FIRE). The lower efficiency emerges due to the coexistence of atoms and nodes with distinct length scales within the multiscale framework, yet the current EMAs fail to differentiate between them. It results in disparate convergence rates across different scales, which undermines both computational accuracy and efficiency. To address the issue, a multiscale FIRE algorithm which updates positions of atoms and nodes synchronously by employing appropriate effective mass is proposed. The optimal effective mass is determined by synchronizing the vibration of harmonic oscillators across different scales. By employing the multiscale FIRE algorithm, the computational efficiency increased by 24.4 and 23.7 times compared to the CG and LBFGS algorithms when used for multiscale nanoindentation simulations. These findings and the proposed algorithm provide valuable insights for structural relaxations of multiscale physical problems and are promising to further improve the computational accuracy and efficiency of MMMs.

Periodic boundary condition along a dislocation line is commonly used in computing activation barriers or formation energies of kink pair of screw dislocation in BCC metals. Although the effect of periodic image interactions on the computation results is obvious, there had been no comprehensive analysis on such effect. In this work, we quantify it through combined nudged elastic band (NEB) simulations and theoretical analysis based on dislocation mechanics. The NEB calculation result demonstrates a non-negligible size dependence on the activation barrier at zero and low stresses. The theoretical analysis offers a practical approach to quantify such size effect without the need of time-consuming NEB simulations. Notably, a simple relationship between kink activation barrier and dislocation line length is derived at zero stress, offering a new approach to compute kink pair formation energy based on NEB simulation results.