Computational insights into the corrosion behavior of NbTaMoW and NbTaMoWV high-entropy alloys in molten fluoride salts

Abstract

Molten salt reactors (MSRs) expose structural materials to harsh conditions, such as elevated temperatures, corrosive fluoride salts, and substantial neutron irradiation. These factors contribute to intricate degradation processes, including radiation-induced defect development, void swelling, and corrosion. Refractory high-entropy alloys with a body-centered cubic structure provide noteworthy thermal stability and mechanical strength, making them excellent candidates for MSR application. This study explores the corrosion properties of NbTaMoW and NbTaMoWV in FLiBe molten salt via density functional theory and ab initio molecular dynamics simulations. Analyses of electronic structure, including density of states and crystal orbital Hamilton population, shed light on interfacial bonding and charge distribution. NbTaMoW shows minimal d-band shifts and weak fluorine interaction, indicating enhanced oxidation resistance. Adding vanadium to form NbTaMoWV further diminishes oxidative vulnerability and stabilizes the electronic structure at the salt interface, suggesting superior corrosion resistance in molten salt conditions.

---

1 Introduction

1.1 High-temperature nuclear energy and the role of molten salts

Nuclear reactors represent the only carbon-neutral energy source capable of providing the multi-gigawatt power generation necessary to satisfy the escalating global electricity demands presented by AI data centers and communities, whilst concurrently facilitating hydrogen production, desalination, and industrial manufacturing.^1,2 Molten salt reactors (MSRs), an ambitious and promising variant of generation-IV nuclear reactors, are distinguished by their inherent passive safety, enhanced thermal efficiency, and effective fuel utilization, employing molten fluoride or chloride salts as a coolant and/or a fuel solvent.³ In contrast to water-cooled reactors, MSRs operate at temperatures ranging from 600 °C to 900 °C and at ambient pressure, owing to the extensive liquid phase characteristics of molten salts. Fluoride salts such as 2LiF–BeF₂ (FLiBe) and LiF–NaF–KF (FLiNaK), offer high heat capacity, thermal conductivity, and stability, making them ideal for use in thermal neutron MSRs.^4,5

1.2 Challenges in molten salt corrosion

Notwithstanding their benefits, MSRs create extremely challenging conditions for structural materials, characterized by elevated temperatures, corrosive molten salts, and substantial neutron irradiation. These parameters provoke intricate degradation processes, including the evolution of radiation-induced defects, void swelling, and corrosion-accelerated damage at grain boundaries.⁶ Alloys that exhibit favorable performance in aqueous settings generally fail in molten salt environments due to the instability of the oxide passivation layer.^7,8 Therefore, the identification of corrosion-resistant materials for fluoride salts is vital for the development of next-generation reactors. Stainless steels and nickel-based alloys are susceptible to degradation through thermodynamic instability, impurity-driven corrosion, selective leaching, and grain boundary attack.^9–13 The Gibbs free energy associated with fluoride formation facilitates metal dissolution, predominantly favoring aluminum, titanium, manganese, and chromium fluorides,¹¹ which become soluble in molten salts, particularly in the absence of an equilibrium.^10,12 The presence of oxygen or moisture intensifies system destabilization, thus complicating efforts in corrosion control. In nickel- and iron-based alloys, a major problem is the loss of chromium, which often starts at the grain boundaries because fluorine atoms gather on the surface.^7,14,15 While HASTELLOY® N performs more effectively than stainless steels at 600–700 °C, it does not possess sufficient mechanical strength at elevated temperatures.¹⁰

1.3 High-entropy alloys – a promising alternative

Alloys suitable for use in molten salt reactors (MSRs) must demonstrate capabilities such as high-temperature durability, resistance to corrosion, tolerance to irradiation, and stable microstructures. High-entropy alloys (HEAs) offer fresh possibilities for cutting-edge nuclear uses by integrating multiple elements in nearly equal proportions, which boosts configurational entropy and prevents the formation of intermetallic phases.¹⁶ Single-phase solid solutions, when stabilized, improve strength, ductility, and resistance to corrosion through lattice distortion, slow diffusion, and compositional customization.^17–19 Materials with body-centered cubic (BCC) structures generally show superior resistance to volumetric swelling compared to face-centered cubic (FCC) varieties, although the specific composition is a key factor. Refractory metals such as Mo, W, and Nb are well known for their outstanding corrosion resistance and reliable operation at high temperatures.^20–22 Refractory metal high-entropy alloys (RHEAs), like NbTaMoW and NbTaMoWV, preserve a single-phase BCC structure at temperatures as high as 1400 °C, showing impressive compressive strength and transitions from ductile to brittle states above room temperature.²³ These characteristics make them ideal for reactor-level thermal fluctuations and structural roles demanding both mechanical robustness and oxidation resistance.

Preliminary investigations have demonstrated that HEAs exhibit enhanced corrosion resistance in aqueous environments when compared to conventional steels.²⁴ A comprehensive review of their aqueous corrosion behavior was recently published by Qiu et al.²⁵ Among the pioneering experimental studies, Chen et al.²⁴ investigated the corrosion resistance of Cu_0.5NiAlCoCrFeSi in H₂SO₄ and NaCl, revealing lower corrosion rates but also a reduced pitting resistance relative to Type-304 stainless steel (304S). The focus has since shifted towards the influence of specific alloying elements. For example, Al contributes to reduced density and increased strength, although it may compromise corrosion resistance due to selective dissolution.²⁶ Ti enhances passivity through the formation of TiO₂ films but may induce heterogeneous microstructures due to intermetallics.²⁷ Chromium facilitates passivation in HEAs by forming Cr₂O₃ films when present in concentrations exceeding 12%, thereby improving corrosion resistance.^7,8,28 Certain chromium-rich alloys, such as Co_1.5CrFeNi_1.5Ti_0.1, demonstrate significant pitting resistance in NaCl solutions,²⁹ although the underlying mechanisms are not yet fully understood. Mo improves both passivity and pit repassivation, but may lead to the formation of sigma phases.³⁰ Ni augments the stability of the FCC structure alongside corrosion resistance but might adversely affect performance due to lattice distortion.³¹

Recent studies have enhanced the understanding of non-aqueous settings. In the work of Elbakhshwan et al.,³² it was found that CrMnFeN high-entropy alloys (HEAs) exposed to molten FLiBe for 1000 hours at 700 °C suffered more mass loss than 316H stainless steel, due to Mn dissolving and the presence of unstable protective oxides similar to Cr₂O₃. This suggest that adjusting the composition could curb degradation. Gwalani et al.³³ utilized atom probe tomography (APT), transmission electron microscopy (TEM), and X-ray absorption near edge structure (XANES) techniques to show that oxidation in CoCrFeNiMn HEAs shifts from atomic-level control to compositional inversion during redox changes, highlighting the complexity of oxide development in harsh conditions. These observations underscore that the corrosion resistance of HEAs fundamentally depends on the choice and proportions of their constituent elements.

1.4 Towards mechanistic and computational insights into HEA corrosion in molten salts

Molten salts present distinct corrosion challenges as compared to aqueous environments, leading to the degradation of conventional alloys through mechanisms such as halide attack, impurity reactions, and selective leaching.^7,8 The utilization of HASTELLOY® N during the 1950s in the Aircraft Reactor Experiment underlines persistent issues associated with corrosion resistance.³⁴ Research indicates that traditional protective oxides become destabilized in environments containing fluoride salts.³² Despite growing interest in HEAs, existing studies predominantly focus on aqueous solutions, which differ from molten salts in terms of electrostatic interactions, ionic strength, and chemical speciation. Techniques such as potentiodynamic polarization (PDP) and electrochemical impedance spectroscopy (EIS) provide only macroscopic insights, leaving the mechanisms at the atomic scale unresolved.^35,36

Recent computational investigations are elucidating the intricate, non-equilibrium characteristics of corrosion in molten salts. The ab initio molecular dynamics simulations conducted by Schneider et al. demonstrated that water impurities in molten NaF and NaCl predominantly remain undissociated, initiating localized oxidation through redox reactions at the metal–salt interface.³⁷ These phenomena are crucial for FLiBe salts, in which trace moisture is inevitable. At a mesoscopic scale, Arkoub et al. employed ReaxFF³⁸ molecular dynamics simulations to examine NiCr alloys in FLiNaK, effectively capturing the kinetics of dissolution and the formation of double-layer structures.³⁹ Xiao et al. utilized kinetic models to elucidate the influence of temperature, voltage, and local chemistry on the processes of oxidation, dissolution, and reprecipitation of Ni–Cr alloys.⁴⁰ Machine learning (ML) frameworks integrating surface energy, phase stability, and Pilling–Bedworth ratios have facilitated the high-throughput screening of corrosion-resistant materials. Random forest models, trained on density functional theory and experimental data, have identified robust HEAs within the AlCrFeCoNi system.⁴¹ These advancements underscore the imperative for ongoing theoretical studies, which integrate atomic-scale simulations, kinetic modeling, and machine learning-driven design, to predict corrosion behaviors in HEAs within complex molten salt environments.

1.5 Objective of this study

In order to address the substantial gap in knowledge concerning the corrosion of high-entropy alloys (HEAs) in molten salt environments, we utilize first-principles methodologies, specifically density functional theory (DFT)⁴² and Ab Initio Molecular Dynamics (AIMD),⁴³ to explore the behavior of two BCC refractory HEAs, NbTaMoW and NbTaMoWV, in a 66–34 mol% LiF–BeF₂ (FLiBe) system. By scrutinizing electronic structure features, such as the density of states (DOS), Crystal Orbital Hamilton Population (COHP) and Electron Localization Function (ELF), we investigate metal–salt interfacial bonding, redox activity, and structural destabilization. This study represents one of the earliest ab initio level investigations into the corrosion behavior of HEAs in molten salt environments, offering mechanistic insight into the initial stages of material degradation.

2 Methods

2.1 Computational framework

Our workflow for designing the interface between the two BCC refractory HEAs, NbTaMoW and NbTaMoWV, and the FLiBe salt is presented in Fig. 1. The HEA and salt subsystems were independently equilibrated using their respective external force fields, 2b+3b tabGAP^44,45 and moment tensor potentials,^46,47 respectively, implemented within the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS).⁴⁸ Equilibration was performed in the isothermal–isobaric (NPT) ensemble at 800 K for both systems, employing a Nose–Hoover thermostat to facilitate thermalization and allow for volume relaxation under ambient pressure. This equilibration step was essential to obtain realistic atomic configurations and densities, especially for the molten salt environment. As illustrated in Fig. 2, a 100 picoseconds (ps) NPT simulation was carried out on BCC (100) slabs of NbTaMoW and NbTaMoWV using 2b+3b tabGAP machine-learned interatomic potentials to generate thermally equilibrated surface geometries.^44,45 In parallel, the Packmol package⁴⁹ was used to create the initial randomized configuration of 66–34 mol% LiF–BeF₂ with 56 F, 28 Li and 14 Be atoms based on an experimental density of 2.156 g cm⁻³ at 800 K.⁵⁰ A separate NPT run was conducted on the salt system using moment tensor machine-learned interatomic potentials,^46,47 matching the in-plane (xy) area of the corresponding HEA slab to ensure interfacial commensurability for subsequent simulations. The HEAs utilized are BCC (100) slabs of (4 × 4 × 5) NbTaMoW and NbTaMoWV, containing a total of 80 atoms (Fig. 2). A vacuum spacing of 20 Å was implemented perpendicular to the surface in the z-direction for both isolated slab and slab–salt interface configurations to mitigate spurious image interactions under periodic boundary conditions. The bottom two layers of each slab were immobilized, whereas the topmost three layers were permitted to move freely.

Fig. 1 Workflow for constructing the Refractory High-Entropy Alloy (RHEA)–FLiBe interface. (Top left) Initial configuration of FLiBe at 0 K. (Top right) Thermally equilibrated FLiBe at 800 K under NPT conditions. (Bottom left) Initial BCC (100)-NbTaMoW slab at 0 K. (Bottom middle) Thermally relaxed BCC (100)-NbTaMoW at 800 K under NPT ensemble. (Bottom right) Construction of the interface between the equilibrated BCC slab and molten salt, under NVT conditions at 800 K.

Fig. 2 Bird's eye view of surface atoms (top) and side view (bottom) of (100)-NbTaMoW and (100)-NbTaMoWV slabs. An NPT simulation was conducted on both (100)-RHEA slabs for 100 ps at 800 K to equilibrate the slab volume and atomic configuration at the given temperature. Only the surface atoms are depicted in the bird's eye view (top) with the atoms in the bottom layers filtered out from the image.

For subsequent first-principles simulations, the equilibrated structures were introduced into Vienna Ab initio Simulation Package (VASP)⁴³ for further refinement through Born–Oppenheimer ab initio molecular dynamics (AIMD). These simulations were executed in a canonical NVT ensemble employing the Nosé–Hoover thermostat at a temperature of 800 K, utilizing a timestep of 1 femtosecond. The electronic structure calculations were executed using spin-polarized density functional theory (DFT) in conjunction with the projector augmented wave (PAW) method. The exchange–correlation functional was addressed within the framework of the generalized gradient approximation (GGA), specifically employing the Perdew–Burke–Ernzerhof (PBE) formulation. For both metal and metal–salt systems, a Γ-centered 1 × 1 × 1 k-point mesh was utilized. A plane-wave energy cutoff of 400 eV was maintained, alongside an electronic convergence threshold of 10⁻⁴ eV. Gaussian smearing with a width of 0.2 eV was implemented for partial occupancies. Van der Waals interactions, which is crucial for accurately modeling non-covalent interactions within the salt and at interface regions, were accounted for using Grimme's D3 dispersion correction.⁵¹ The simulation cells comprised of isolated HEA slabs (NbTaMoW or NbTaMoWV), including HEA slabs interfaced with FLiBe. Each system underwent equilibration at 800 K for no less than 10 ps, followed by production runs of up to 30 ps for the purposes of statistical averaging.

To validate the reliability of our DFT framework and the accuracy of the interatomic potentials used, we performed a series of equation of state (EOS) calculations for each pure BCC refractory metal (Nb, Ta, Mo, W, V) using both DFT and the 2b+3b tabGAP machine-learned potential.⁴⁴ These calculations were conducted within the Atomic Simulation Environment (ASE) framework, using the VASP and LAMMPS calculators, respectively. Murnaghan EOS models⁵² were constructed from static energy calculations at fixed volumes, and key material properties—including lattice parameter and bulk modulus—were extracted and compared against experimental values.⁵³ The DFT predictions showed excellent agreement with experimental lattice constants (within 0.3% error, except for V), and reasonable accuracy for bulk moduli which agree with the results of Fernández-Pello et al.⁵⁴ The tabGAP potential also reproduced lattice constants with good fidelity, supporting its transferability to pure-element systems. Full details and comparative results, including EOS fits and property discrepancies, are provided in the supplementary information (SI) (Fig. S3, S4 and Table S3).

2.2 Electronic structure and bonding analysis

The electronic structure and bonding characteristics were investigated utilizing Density of States (DOS), Crystal Orbital Hamilton Population (COHP) analyses⁵⁵ and the partitioning of the Electron Localization Function (ELF) and subsequent integration of electron density.⁵⁶ DOS analysis provides insights into the distribution of electronic states across energy levels, facilitating the identification of material properties, bonding characteristics and the contribution of specific atomic orbitals to the electronic structure. To examine covalent bonding in real space, COHP analysis was employed using the LOBSTER code,⁵⁷ which reconstructs localized orbital information from plane-wave DFT outputs. COHP divides the band-structure energy into bonding, nonbonding, and antibonding contributions between atom pairs, which are projected as energy-resolved curves. Negative COHP values indicate bonding states, while positive values denote antibonding interactions. This methodology establishes a direct correlation between the electronic structure and chemical bonding, aiding in the identification of key bonding interactions, bond strengths, and their energetic contributions to phase stability. For the evaluation of electron pairing and bonding localization, the Electron Localization Function (ELF) was computed and analyzed using Critic2,⁵⁶ a real-space post-processing tool. ELF partitioning was used to reveal the charge distribution across the metal–salt interface and to assess the impact of salt contaminants on local electronic environments.

3 Results

3.1 HEA subsystem equilibration

Fig. 2 illustrates that both (100)-NbTaMoW and NbTaMoWV slabs are no longer a single-phase BCC crystalline after structure relaxation at 800 K. Many atoms have moved away from their perfectly symmetrical BCC lattice arrangement at 0 K to optimize the energy of the structure at 800 K. In particular, a common neighbor analysis of the equilibrated (100)-NbTaMoW and NbTaMoWV slabs highlights an ≈50% loss of the BCC crystal structure at 800 K (Fig. S7). This behavior is attributed to the formation of point defects, such as metallic vacancies and self-interstitials, and the enhanced motion of dislocations at elevated temperatures, which become amplified in smaller alloy systems. Variation in the energy landscape is more likely to be observed in complex, multi-component alloys whose constituent elements have vastly differing Peierls barriers.⁵⁸ Hence, the movement of dislocations in NbTaMoW and NbTaMoWV becomes unpredictable, and the refractory HEAs adopt an asymmetric (or lopsided) configuration due to kinks in the lattice as seen in Fig. 2.

3.2 Density of states

In Fig. 3 and 4, we present the DOS of NbTaMoW and NbTaMoWV slabs, both in vacuum and in contact with molten FLiBe. These configurations are analyzed at 40 ps of simulation time, respectively. Two simulations were conducted in parallel, with one including the salt overlay represented by the red line, and the other without the salt represented by the blue line. To ensure a direct comparison between the alloy surfaces devoid of salt and those covered with salt, the energy axis remains unadjusted relative to the Fermi level (E_F). Fig. 5 compares the NbTaMoW–FLiBe and NbTaMoWV–FLiBe configurations after 40 ps, highlighting the similarities and differences in the electronic stability of the two RHEA compositions in the presence of salt. The salt–alloy interfaces without and with vanadium (V) are represented by the green and purple line, respectively. The DOS in proximity to E_F is particularly significant due to its profound influence on bonding characteristics and material properties.

Fig. 3 Total and partial density of states of the BCC (100)-NbTaMoW surface, in contact with a molten FLiBe salt after 40 ps at 800 K. The blue curves represent the DOS of the clean refractory high-entropy alloy (RHEA) surface without the salt layer, while the red curves show the DOS with the FLiBe layer present. The dotted blue and red vertical lines represent the Fermi level of the clean RHEA surface and the RHEA with the added salt, respectively.

Fig. 4 Total and partial density of states of the BCC (100)-NbTaMoWV surface, in contact with a molten FLiBe salt after 40 ps at 800 K. The blue curves represent the DOS of the clean refractory high-entropy alloy (RHEA) surface without the salt layer, while the red curves show the DOS with the FLiBe layer present The dotted blue and red vertical lines represent the Fermi level of the clean RHEA surface and the RHEA with the added salt, respectively.

Fig. 5 Total and partial density of states of BCC (100)-NbTaMoW and NbTaMoWV surfaces, in contact with a molten FLiBe salt after 40 ps at 800 K. The green curves represent the DOS of the (100)-NbTaMoW–FLiBe interface, while the purple curves represent the DOS of the (100)-NbTaMoWV–FLiBe interface. The dotted green and purple vertical lines represent the Fermi level of the NbTaMoW and NbTaMoWV alloy with the added salt, respectively.

As anticipated for metallic systems, both the total and partial DOS of NbTaMoW and NbTaMoWV reveal substantial overlap between valence and conduction bands, thereby substantiating the metallic nature of these refractory high-entropy alloys (RHEAs). This overlap is evident in Fig. 3 and 4. In both systems, irrespective of the presence of a FLiBe layer (red line), the total DOS manifests a bimodal distribution close to the Fermi level, with the Fermi level located within a “pseudo-band gap”. This feature is characteristic of the stable BCC structures in RHEAs, as previously documented.^59,60 Despite the formation of defects in the RHEAs at 800 K, there is still some BCC crystalline character remaining (see Fig. S7) which is confirmed by the “pseudo-band gap” in the DOS. Nevertheless, the introduction of the FLiBe salt layer (red line) markedly broadens this pseudo-band gap in the NbTaMoW system relative to NbTaMoWV which becomes apparent in Fig. 5. This phenomenon is consistent with the findings of Hu et al.⁶¹ and implies an enhancement in covalent character among Nb, Ta, Mo, and W atoms in the presence of the salt. Conversely, the incorporation of vanadium (V) appears to enhance metallic bonding within the alloy, thereby diminishing the pseudo-band gap and elucidating its role in modulating the electronic structure.

For both HEA–FLiBe interfaces, the primary constituents of the lower valence band are the 2p-orbitals of F, whereas the upper valence band and lower conduction band are predominantly composed of metal d-orbitals. This electronic configuration suggests an enhanced covalent character associated with the 2p-orbitals of F and metallic bonding due to the delocalized nature of the d electrons. Furthermore, the detected overlap between metal s and d orbital contributions within the valence band signifies s–d hybridization, which is a hallmark of metallic cohesion in transition metals. The minimal overlap between the 2p of F and metal d states also suggests weak orbital hybridization between F and the metal sublattice, indicating modest chemical interaction at the interface. Significantly, the partial d-orbital DOS for Nb and Ta (see Fig. 3 Nb, Ta) manifests a left-skewed distribution with less pronounced bimodal features than those observed in Mo and W (Fig. 3 Mo, W). This asymmetry indicates that Nb and Ta exhibit a higher d-orbital DOS in the conduction band compared to the valence band, suggesting a lower electronic stability. Conversely, Mo and W exhibit a nearly symmetric bimodal d-band structure, typically associated with greater electronic stability and a reduced propensity for oxidation during the degradation processes of HEAs. This enhanced stability found in Mo and W arises from their half-filled d-orbitals, which elevate ionization energies and diminish their reactivity in corrosive environments.

As seen in Fig. 3 and 4, upon exposure to FLiBe, both HEAs exhibit a shift in their occupied electronic states toward higher energy levels, resulting in a progressive increase in E_F. This evolution is captured in DOS at 20 ps and 30 ps, shown in Fig. S1 and S2 in the SI. At 40 ps, the separation between the 2p orbitals of F and the d orbitals of the metal has increased, due to the upward shift of the metal d-states and a downward shift of the F 2p-states – indicative of progressive oxidation by fluorine. Among the constituent elements, Mo d-electrons show the smallest energy shift, suggesting lower susceptibility to oxidation, while Ta d-electrons exhibit the largest shift. Notably, the NbTaMoW–FLiBe interface shows a greater increase in high-energy electron occupation compared to the NbTaMoWV–FLiBe interface (Fig. 5), indicating that the addition of V enhances corrosion resistance. This is consistent with the stronger metal–metal bonding introduced by vanadium in the NbTaMoWV HEAs.

3.3 Crystal orbital Hamilton population analysis

To gain quantitative insight into the bonding interactions between fluorine species and refractory metal atoms at the salt–HEAs interface, Crystal Orbital Hamilton Population (COHP) analysis was performed on both NbTaMoW–FLiBe and NbTaMoWV–FLiBe systems. The analysis focused on time snapshots taken at 10, 20, 30, and 40 ps during molecular dynamics simulations conducted at 800 K. The evolution of bond lengths and corresponding integrated COHP (ICOHP) values for selected representative metal–fluorine (M–F) pairs with distances below 3 Å are presented in Fig. 6 and 7. The electronic structure is illustrated by the COHP curves shown in Fig. 8 and further detailed in Tables S1 and S2, as well as Fig. S5 and S6 in the SI.

Fig. 6 Bond length and integrated crystal orbital Hamilton population (ICOHP) of metal fluorides at the (100)-NbTaMoW–FLiBe interface at 800 K. The cutoff distance for the atom pair is 3 Å. If the distance between the atom pair exceeds 3 Å, the interatomic interaction is considered insignificant and the bond length and ICOHP are not reported.

Fig. 7 Bond length and integrated crystal orbital Hamilton population (ICOHP) of metal fluorides at the (100)-NbTaMoWV–FLiBe interface at 800 K. The cutoff distance for the atom pair is 3 Å. If the distance between the atom pair exceeds 3 Å, the interatomic interaction is considered insignificant and the bond length and ICOHP are not reported.

Fig. 8 Crystal orbital Hamilton population (COHP) analysis of the BCC (100) surfaces of NbTaMoW (top) and NbTaMoWV (bottom) in contact with molten FLiBe, evaluated at 800 K after 40 ps of simulation time. The black dotted line represents the Fermi level which is set to 0 eV.

Across both systems and all time steps (Fig. 6 and 7), no persistent covalent M–F bonds were observed. Specifically, bond lengths frequently exceed 2.0 Å and ICOHP values often approach zero, indicating minimal orbital overlap and hence negligible chemical bonding. This result strongly suggests that the corrosive interaction between molten FLiBe and the HEAs surfaces is minimal under these conditions, with no evidence of sustained M–F covalent bonding between 10 and 40 ps at 800 K. These findings reinforce the chemical inertness of the NbTaMoW and NbTaMoWV surfaces in contact with pure molten FLiBe salt.

Nonetheless, transient bonding events are observed at specific time frames. For the NbTaMoW–FLiBe interface at 20 ps (Fig. 6), the strongest M–F interactions were observed between Ta20–F3 and Ta8–F46, reflected in shorter bond lengths and moderately negative ICOHP values (e.g., −1.5 to −2.0 eV). At 30 ps, the interaction between Nb20–F14 became dominant, characterized by a bond length of 1.781 Å and a significantly negative ICOHP value of −2.345 eV – the most negative of all pairs examined – indicating that Nb exhibits the strongest affinity for fluorine among the constituent metals. At 40 ps, a relatively strong interaction was also observed between W2–F6, although with a less negative ICOHP compared to Nb–F, suggesting a weaker but more stable bonding tendency.

Significantly, it is evident that M–F interactions at distances less than 2 Å frequently occur in the NbTaMoWV–FLiBe interface (Fig. 7) as compared to NbTaMoW–FLiBe (Fig. 6), although these additional interactions do not necessarily confer enhanced strength. Notably, V–F bonds frequently manifest at the NbTaMoWV interface; however, they generally display relatively weak interaction strengths, as evidenced by modest ICOHP values. For instance, at 10 ps, the V3–F16 and Ta4–F47 pairs form bonds with distances shorter than 2 Å. The V3–F16 bond persists throughout the simulation up to 30 ps, maintaining a consistent bond length between 1.90 and 2.00 Å, which suggests spatial stability. Though less enduring, Nb2–F42 also retains bonding from 10 to 30 ps, achieving a peak ICOHP of −1.251 eV at 20 ps. Among the more robust transient interactions, Ta–F bonds are particularly noteworthy: F47–Ta4 (10–30 ps) and F24–Ta1 (20 ps) exhibit significant overlap. However, the strengths of Ta–F bonds experience considerable fluctuations over time, indicating a dynamic process of bond formation and dissolution that reflects interfacial mobility rather than enduring chemisorption.

In Fig. 8, the electronic structure illustrated by the COHP curves is presented for selected metal–fluorine bonds at the BCC (100) surfaces of NbTaMoW (top, 40 ps) and NbTaMoWV (bottom, 40 ps) in contact with molten FLiBe at 800 K. In the NbTaMoW interface (40 ps), three significant M–F interactions are discerned: Ta20–F3, Nb20–F14, and W2–F6. Among these, W2–F6 demonstrates the most prominent bonding peak centered around −7 eV, indicative of a relatively strong covalent character. The COHP characteristics for Ta20–F3 and Nb20–F14 are also of considerable magnitude, implying meaningful, albeit slightly weaker bonding. These results underscore Nb and Ta as crucial bonding centers within the pristine alloy interface and are consistent with prior ICOHP and ELF data, which indicate their reactivity toward fluorine. Conversely, the NbTaMoWV–FLiBe interface at 40 ps exhibits a more intricate bonding environment, characterized by a broader array of M–F interactions involving V, Ta, Nb, and Mo. Notably, Ta4–F47 and Ta1–F24 manifest robust bonding features, including pronounced −COHP peaks below the Fermi level, aligning with transient yet significant metal–fluorine interactions. V3–F16 and V15–F52 also illustrate moderate bonding contributions, affirming the role of vanadium in modulating interfacial reactivity, albeit with less pronounced COHP intensities. Intriguingly, Mo3–F33 exhibits a sharp bonding peak around −6 eV, whereas Nb2–F42 demonstrates weak contribution. Overall, the presence of vanadium augments the diversity of M–F interactions at the interface; however, the bonding is more dispersed and generally less robust compared to the concentrated and stronger interactions within the NbTaMoW system. This observation corroborates the interpretation that while V-rich interfaces may facilitate a higher frequency of bonds with fluorine, they do not necessarily promote stronger or more stable interactions. These findings provide insights into the bonding dynamics and potential corrosion resistance of HEAs systems.

In summary, analyses of COHP and ICOHP elucidate the following findings: (i) stable covalent metal–fluoride bonds are absent throughout the 10–40 ps simulation interval at 800 K, (ii) among the transient bonds, Nb–F interactions emerge as the most robust, whereas V–F bonds occur with greater frequency but exhibit reduced strength, (iii) the NbTaMoWV alloy system manifests a higher incidence of weak M–F contacts, which suggests an increase in fluorine coordination, albeit without corresponding enhancement in interfacial bonding strength, (iv) Ta–F bonds demonstrate notable dynamism, intermittently forming and disbanding, indicative of a fluid and weakly interactive interface. These ephemeral yet measurable interactions suggest that, although the majority of M–F contacts are weak or non-bonding, particular local geometries – especially those involving Nb and W – have the potential to temporarily facilitate bonding. The lack of enduring covalent M–F bonds across all observed time frames and alloy compositions highlights the limited corrosive reactivity of FLiBe with RHEA surfaces at elevated temperatures, in the absence of impurities such as moisture or oxygen.

3.4 Electron localization function

In order to elucidate the interfacial redox behavior at elevated temperatures, Electron Localization Function (ELF) analysis was conducted for both NbTaMoW–FLiBe and NbTaMoWV–FLiBe interfaces over a simulation duration ranging from 10 to 40 ps, with the number of valence electrons calibrated against the valence electrons in the initial structure (at 0 ps). The corresponding change in valence electrons with time for NbTaMoW–FLiBe and NbTaMoWV–FLiBe are detailed in Fig. 9 and 10. Within this context, a reduction in electron density signifies oxidation, which is commonly observed on transition metals, whereas an increase indicates reduction, predominantly occurring on electronegative fluorine atoms. As reported in Section 3.3, at the NbTaMoW interface, notwithstanding that the Nb20–F14 interaction demonstrates the lowest ICOHP (Fig. 6) value at 30 ps, the Ta20–F3 bond exhibits more pronounced polarization. Nonetheless, both bonds undergo destabilization over time, culminating in the dissociation of Ta20–F3 by 30 ps and Nb20–F14 by 40 ps. This observation is supported by the evolution of ELF, which signifies considerable electron density restoration on the metal sites and depletion on the fluorine atoms, with Nb20 and F14 reverting to their initial charge states, indicating complete bond reversal. Importantly, despite the robust initial bonding, Nb20–F14 does not maintain a low reduction potential over an extended period, highlighting that transient bonding strength does not inherently equate to enduring redox stability.

Fig. 9 Change in valence electrons of surface atoms at the interface between (100)-NbTaMoW and FLiBe at 800 K. The number of valence electrons for a given atom is normalized to the valence electrons at 0 ps.

Fig. 10 Change in valence electrons of surface atoms at the interface between (100)-NbTaMoWV and FLiBe at 800 K. The change in valence electrons for atoms participating in V–F bonds is at the top and the atoms involved in other metal fluoride bonds are at the bottom. The number of valence electrons in a given atom is normalized to the valence electrons at 0 ps.

At the NbTaMoWV interface (Fig. 10), similar trends are observed with notable differences due to the presence of vanadium. For instance, the V2–F3 bond undergoes rapid dissociation at 10 ps, as indicated by a sharp increase in valence electrons in V2 +0.38 (check V2 in ELF plot), confirming oxidation-induced bond rupture. While V3–F16 and V15–F52 appear more geometrically stable, ELF data show progressive oxidation of V3 and V15, suggesting that V participates in dynamic but weak bonding interactions with fluorine. Additionally, Ta1 at the NbTaMoW interface shows the most pronounced charge depletion at 20 ps, aligning with the peak in Ta1–F24 bond strength (Fig. 7), and indicating a strongly polarized and chemically active site. Across both interfaces, tantalum consistently exhibits significant electron loss, pointing to a higher oxidation susceptibility relative to other constituent metals. This trend is consistent with the elevated electronic energy near the Fermi level observed in the Ta d-orbital projected DOS, which reflects a higher tendency for Ta to donate electrons to electronegative species such as fluorine.

3.5 Effect of RHEAs on salt chemistry

The radial distribution function was calculated for FLiBe salt after contact with both (100)-NbTaMoW and -NbTaMoWV slabs at 800 K to deduce the effect of the RHEAs on the salt chemistry (Fig. 11). The structure of the salt remains largely unchanged after exposure to the alloy surfaces for 40 ps at 800 K, which is expected provided the short simulation time. The salt maintains its liquid phase, as depicted by the short-range ordering in Fig. 11. Both salt–alloy interfaces exhibit a small decrease in the Be–F bonding strength after 40 ps. Furthermore, there is an increase in the F–F bonding strength in NbTaMoW–FLiBe and a decrease in the F–F bonding strength in NbTaMoWV–FLiBe. The shift in the Li–F bonding strength is insignificant at both interfaces.

Fig. 11 Radial distribution function of FLiBe salt after contact with (100)-NbTaMoW (top) and -NbTaMoWV (bottom) at 800 K for 40 ps.

4 Conclusions

In this study, the corrosion behavior of NbTaMoW and NbTaMoWV HEAs in molten salt environments was examined using Density Functional Theory (DFT) and Ab Initio Molecular Dynamics (AIMD) simulations. Electronic structure analyses, comprising Density of States (DOS), Crystal Orbital Hamilton Population (COHP), and Electron Localization Function (ELF) partitioning, were employed to reveal the bonding characteristics and charge distribution at the interface. Our key findings are summarized as follows: (i) Element-specific oxidation behavior: Ta is the most susceptible to oxidation, as evidenced by significant shifts of d-orbitals to higher energies, unstable Ta–F bonds (as shown by COHP), and substantial electron transfer to fluorine (confirmed by ELF analysis). Mo, in contrast, displays the smallest shift in d-states and minimal interaction with fluorine, suggesting greater resistance to oxidation. Nb forms the strongest interactions with fluorine F in NbTaMoW, as indicated by the short bond length and highly negative ICOHP. (ii) Effect of vanadium addition: V forms weak bonds with fluorine (ICOHP > −1.0 eV), indicating minimal redox activity and low corrosion reactivity. The integration of vanadium into the NbTaMoW alloy (i.e., forming NbTaMoWV) results in reduced electronic excitation and weaker metal–fluorine interactions, suggesting enhanced structural and chemical stability. (iii) Interface stability: the NbTaMoW interface shows stronger Nb–F interactions and greater d-electron energy shifts, indicating a higher tendency towards oxidation compared to the NbTaMoWV interface. This observation reinforces the conclusion that vanadium enhances corrosion resistance by strengthening metal–metal bonding and reducing the availability of reactive sites. (iv) Salt chemistry: the structure of the FLiBe salt remains largely unaffected in the presence of NbTaMoW and NbTaMoWV after 40 ps at 800 K. Overall, the incorporation of vanadium into NbTaMoW-based HEAs appears to mitigate corrosion by reducing oxidation susceptibility and stabilizing the alloy's electronic structure at the molten salt interface. These insights highlight the potential of tailored compositional adjustments in HEAs for improved performance in extreme environments.

Author contributions

A. B. K. and C. G. T. F designed the computational framework. A. B. K. implemented the framework, performed all calculations, and analyzed the data. A. B. K. drafted the manuscript with input from all authors. A. V. R. V. assisted in writing – review and editing. C. G. T. F. conceived the study and supervised the overall direction and planning.

Conflicts of interest

There are no conflicts of interest to declare.

Data availability

Data supporting this article have been included as part of the supplementary information (SI). All MD simulations and DFT calculations data are available at https://doi.org/10.5281/zenodo.15699534. Supplementary information: Fig. S1 and S2: DOS of RHEAs with FLiBe salt after 20 and 30 ps. Tables S1 and S2: tabulated bond lengths and ICOHP values for NbTaMoW and NbTaMoWV with FLiBe salt. Fig. S3 and S4: EOS for refractory metal unit cells using DFT with PBE and tabGAP potentials. Fig. S5 and S6: COHP analysis of NbTaMoW and NbTaMoWV with FLiBe salt throughout the simulation. Fig. S7: common neighbor analysis of RHEAs after NPT simulation with classical MD. Table S3: tabulated EOS data for refractory metal unit cells using DFT with PBE and tabGAP potentials. See DOI: https://doi.org/10.1039/d5fd00102a.

Acknowledgements

This work was supported by NSERC-CNSC Small Modular Reactors Research Grant Initiative. This work was made possible by the facilities of the Shared Hierarchical Academic Research Computing Network (SHARCNET: https://www.sharcnet.ca) and Digital Research Alliance of Canada (https://alliancecan.ca/en).

References

1. R. S. El-Emam, A. Constantin, R. Bhattacharyya, H. Ishaq and M. E. Ricotti, Renewable Sustainable Energy Rev., 2024, 192, 114157.
2. S. Aslam, S. Rani, K. Lal, M. Fatima, T. Hardwick, B. Shirinfar and N. Ahmed, Green Chem., 2023, 25, 9543–9573.
3. G.I.F Portal, Molten Salt Reactors (MSR), https://www.gen-4.org/generation-iv-criteria-and-technologies/molten-salt-reactors-msr.
4. J. Barnes, R. Coutts, T. Horne and J. Thai, PAM Review Energy Science & Technology, 2019, 6, 38–55.
5. D. F. Williams, Assessment of Candidate Molten Salt Coolants for the Advanced High Temperature Reactor (AHTR), 2006.
6. Corrosion in Molten Salt Reactors, https://madcor.neep.wisc.edu/previous-work/corrosion-in-molten-salt-reactors, The UW-Madison MAterials Degradation under COrrosion and Radiation (MADCOR).
7. A. Leong, J. Zhang and S. D. Rountree, High Temp. Corros. Mater., 2023, 99, 375–397.
8. A. E. Danon, O. Muránsky, I. Karatchevtseva, Z. Zhang, Z. J. Li, N. Scales, J. J. Kruzic and L. Edwards, Corros. Sci., 2020, 164, 108306.
9. J. DeVan and R. Evans III, ORNL/TM-328, Oak Ridge National Laboratory, Oak Ridge, TN, 1962.
10. J. DiStefano, J. DeVan, J. Keiser, R. Klueh and W. Eatherly, Materials Considerations for Molten Salt Accelerator-Based Plutonium Conversion Systems, Oak Ridge National Laboratory, technical report, 1995.
11. K. Sridharan and T. R. Allen, Molten Salts Chemistry, Elsevier, 2013, pp. 241–267.
12. L. C. Olson, J. W. Ambrosek, K. Sridharan, M. H. Anderson and T. R. Allen, J. Fluor. Chem., 2009, 130, 67–73.
13. S. Guo, J. Zhang, W. Wu and W. Zhou, Prog. Mater. Sci., 2018, 97, 448–487.
14. L. Olson, K. Sridharan, M. Anderson and T. Allen, Mater. High Temp., 2010, 27, 145–149.
15. C.-L. Ren, H. Han, W.-B. Gong, C.-B. Wang, W. Zhang, C. Cheng, P. Huai and Z.-Y. Zhu, J. Nucl. Mater., 2016, 478, 295–302.
16. J.-W. Yeh, S.-K. Chen, S.-J. Lin, J.-Y. Gan, T.-S. Chin, T.-T. Shun, C.-H. Tsau and S.-Y. Chang, Adv. Eng. Mater., 2004, 6, 299–303.
17. E. P. George, D. Raabe and R. O. Ritchie, Nat. Rev. Mater., 2019, 4, 515–534.
18. W.-L. Hsu, C.-W. Tsai, A.-C. Yeh and J.-W. Yeh, Nat. Rev. Chem., 2024, 8, 471–485.
19. J. Liang, W. Wang, Z. Cao, J. Guo, Z. Sun and Y. Hai, Mater. Corros., 2024, 75, 424–432.
20. A. K. Misra and J. D. Whittenberger, 22nd Intersociety Energy Conversion Engineering Conference, 1987, p. 9226.
21. K. M. Sankar and P. M. Singh, Corros. Sci., 2023, 213, 110977.
22. R. Odette and S. Zinkle, Structural Alloys for Nuclear Energy Applications, Newnes, 2019.
23. O. N. Senkov, G. B. Wilks, J. M. Scott and D. B. Miracle, Intermetallics, 2011, 19, 698–706.
24. Y. Chen, T. Duval, U. Hung, J. Yeh and H. Shih, Corros. Sci., 2005, 47, 2257–2279.
25. Y. Qiu, S. Thomas, M. A. Gibson, H. L. Fraser and N. Birbilis, npj Mater. Degrad., 2017, 1, 15.
26. C. P. Lee, C. C. Chang, Y. Y. Chen, J. W. Yeh and H. C. Shih, Corros. Sci., 2008, 50, 2053–2060.
27. A. Mazzarolo, M. Curioni, A. Vicenzo, P. Skeldon and G. Thompson, Electrochim. Acta, 2012, 75, 288–295.
28. K. Sieradzki and R. Newman, J. Electrochem. Soc., 1986, 133, 1979–1980.
29. Y. L. Chou, Y. C. Wang, J. W. Yeh and H. C. Shih, Corros. Sci., 2010, 52, 3481–3491.
30. Y. L. Chou, J. W. Yeh and H. C. Shih, Corros. Sci., 2010, 52, 2571–2581.
31. X.-W. Qiu and C.-G. Liu, J. Alloys Compd., 2013, 553, 216–220.
32. M. Elbakhshwan, W. Doniger, C. Falconer, M. Moorehead, C. Parkin, C. Zhang, K. Sridharan and A. Couet, Sci. Rep., 2019, 9, 18993.
33. B. Gwalani, A. Martin, E. Kautz, B. Guo, S. Lambeets, M. Olszta, A. K. Battu, A. Malakar, F. Yang and J. Guo, et al. , Nat. Commun., 2024, 15, 5026.
34. S. Delpech, E. Merle-Lucotte, D. Heuer, M. Allibert, V. Ghetta, C. Le-Brun, X. Doligez and G. Picard, J. Fluor. Chem., 2009, 130, 11–17.
35. D. Chen, W. Zhou, Y. Ji and C. Dong, MGE Adv., 2025, 3, 83.
36. S. Li, C. Li and F. Wang, Mater. Today Chem., 2024, 37, 101986.
37. A. Schneider, D. Andersson and Y. Zhang, Commun. Mater., 2024, 5, 91.
38. T. P. Senftle, S. Hong, M. M. Islam, S. B. Kylasa, Y. Zheng, Y. K. Shin, C. Junkermeier, R. Engel-Herbert, M. J. Janik, H. M. Aktulga, T. Verstraelen, A. Grama and A. C. T. van Duin, npj Mater. Degrad., 2016, 2, 15011.
39. H. Arkoub, S. Dwivedi, A. C. T. van Duin and M. Jin, Appl. Surf. Sci., 2024, 655, 159627.
40. P. Xiao, C. A. Orme, S. R. Qiu, T. A. Pham, S. Cho, M. Bagge-Hansen and B. C. Wood, Nat. Commun., 2025, 16, 341.
41. C. Zeng, A. Neils, J. Lesko and N. Post, Comput. Mater. Sci., 2024, 237, 112925.
42. W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, 1133–1138.
43. G. Kresse and J. Furthmüller, Phys. Rev. B:Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
44. J. Byggmästar, K. Nordlund and F. Djurabekova, Phys. Rev. B, 2021, 104, 104101.
45. J. Byggmästar, K. Nordlund and F. Djurabekova, Phys. Rev. Mater., 2022, 6, 083801.
46. A. V. Shapeev, Multiscale Modeling and Simulation, 2016, vol. 14, pp. 1153–1173.
47. S. Attarian, D. Morgan and I. Szlufarska, J. Mol. Liq., 2022, 368, 120803.
48. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, Comput. Phys. Commun., 2022, 271, 108171.
49. L. Martínez, R. Andrade, E. G. Birgin and J. M. Martínez, J. Comput. Chem., 2009, 30, 2157–2164.
50. R. R. Romatoski and L. W. Hu, Ann. Nucl. Energy, 2017, 109, 635–647.
51. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Phys. Chem., 2010, 132, 154104.
52. C. L. Fu and K. M. Ho, Phys. Rev. B:Condens. Matter Mater. Phys., 1983, 28, 5480–5486.
53. WebElements, https://www.webelements.com/.
54. D. Fernández-Pello, J. M. Fernández-Díaz, M. A. Cerdeira, C. González and R. Iglesias, Mater. Today Commun., 2020, 24, 101323.
55. R. Dronskowski and P. E. Bloechl, J. Phys. Chem., 1993, 97, 8617–8624.
56. A. Otero-de-la Roza, E. R. Johnson and V. Luaña, Comput. Phys. Commun., 2014, 185, 1007–1018.
57. S. Maintz, V. L. Deringer, A. L. Tchougréeff and R. Dronskowski, J. Comput. Chem., 2016, 37, 1030–1035.
58. P. Borges, R. Ritchie and M. Asta, Sci. Adv., 2024, 10, eadp7670.
59. Y. J. Hu, G. Zhao, B. Zhang, C. Yang, M. Zhang, Z. K. Liu, X. Qian and L. Qi, Nat. Commun., 2019, 10, 4484.
60. A. P. Sutton, Electronic Structure of Materials, Oxford University Press, 1993.
61. Y. L. Hu, L. H. Bai, Y. G. Tong, D. Y. Deng, X. B. Liang, J. Zhang, Y. J. Li and Y. X. Chen, J. Alloys Compd., 2020, 827, 153963.

---