Complex strengthening mechanisms in the NbMoTaW multi-principal element alloy

Refractory multi-principal element alloys (MPEAs) have exceptional mechanical properties, including high strength-to-weight ratio and fracture toughness, at high temperatures. Here we elucidate the complex interplay between segregation, short-range order, and strengthening in the NbMoTaW MPEA through atomistic simulations with a highly accurate machine learning interatomic potential. In the single crystal MPEA, we find greatly reduced anisotropy in the critically resolved shear stress between screw and edge dislocations compared to the elemental metals. In the polycrystalline MPEA, we demonstrate that thermodynamically driven Nb segregation to the grain boundaries (GBs) and W enrichment within the grains intensifies the observed short-range order (SRO). The increased GB stability due to Nb enrichment reduces the von Mises strain, resulting in higher strength than a random solid solution MPEA. These results highlight the need to simultaneously tune GB composition and bulk SRO to tailor the mechanical properties of MPEAs.

Despite intense research efforts, the fundamental mechanisms behind the remarkable mechanical properties of MPEAs remain under heavy debate. Solid solution strengthening, whereby the existence of multiple elements components of different atomic radii and elastic moduli impede dislocation motion, has been proposed as a key mechanism in both fcc and bcc MPEAs 9,10 . However, it is clear that the microstructure (e.g., nanotwinning), short-range order (SRO), phase transitions, and other effects also play significant roles [11][12][13][14] .
Computational simulations are an important tool to elucidate the fundamental mechanisms behind the observed strengthening in MPEAs. However, due to the high computational cost, investigations of MPEAs using density functional theory (DFT) calculations have been limited to bulk special quasi-random structures (SQSs) [15][16][17] . Atomistic simulations using linear-scaling interatomic potentials (IAPs) can potentially access more complex models and longer timescales. However, classical IAPs, such as those based on the embedded atom method, are fitted mainly to elemental properties and generally perform poorly when scaled to multi-component alloys. Furthermore, classical IAPs are typically explicit parameterizations of two-body, three-body, and manybody interactions and hence becomes combinatorially complex for multi-element systems, such as MPEAs 18 . Recently, machine learning of the potential energy surface as a function of local environment descriptors has emerged as a systematic, reproducible, automatable approach to develop IAPs (ML-IAPs) with near-DFT accuracy for elemental as well as multi-component systems [19][20][21][22][23][24][25][26][27] . While a few ML-IAPs have been developed for MPEAs 28 , they have mainly applied to the study of phase stability of the bulk alloy.
In this work, we develop a ML-IAP for the refractory NbMoTaW alloy system using the spectral neighbor analysis potential (SNAP) approach 22 . Using this MPEA SNAP model, we show that the Peierls stress for both screw and edge dislocation in the equiatomic NbMoTaW MPEA are much higher than those for all the individual metals, and edge dislocations become much more important in the MPEA than that in the pure elemental bcc system. From Monte Carlo (MC)/molecular dynamics (MD) simulations, we find strong evidence of Nb segregation to the grain boundaries (GBs) of the NbMoTaW MPEA, which in turn has a substantial effect on the observed SRO. The observed Nb segregation to the GB leads to an enhancement in the strength of the MPEA.

RESULTS
NbMoTaW SNAP model Figure 1 shows the workflow adopted in fitting the quaternary NbMoTaW MPEA SNAP model and methodological details are provided in the "Methods" section. Briefly, the MPEA SNAP model was fitted in three steps, as illustrated by the right panel of Fig. 1 with three optimization units 24,25 from left to right. In the first step, a SNAP model was fitted for each component element, as shown in the left optimization unit of the right panel of Fig. 1. The optimized SNAP model coefficients β are provided in Supplementary Table 1, and the mean absolute error (MAE) in energies and forces are provided in Supplementary Fig. 1. The optimized radius cutoffs R El c for Nb, Mo, Ta, and W are 4.7, 4.6, 4.5, and 4.5 Å, respectively, which are slightly larger than the third nearestneighbor distance for each element. This result is in line with previous models developed for bcc elements 22,24,25,27 . These optimized radius cutoffs were adopted for the MPEA SNAP model fittings in the next two steps. In the second step, the data weights (ω) were fixed according to the number of data points for each data group. A grid search was performed for the atomic weights (S k atom ) by generating a series of SNAP models with different combinations of atomic weights by only running the inner loop in the optimization unit 24,25 , as shown in the middle optimization unit of the right panel of Fig. 1. In the third step (see the right optimization unit of the right panel of Fig. 1), the ten combinations of atomic weights (see Supplementary Table 2) with the best accuracy in energy and force predictions were chosen to conduct a full optimization, including the optimization of the data weight in the outer loop of the optimization unit. The optimized parameters of the best model, i.e., the model with the smallest MAE in energies and forces, are provided in Supplementary Table 3. The training data comprises DFT-computed energies and forces for ground-state structures, strained structures, surface structures, SQSs 29 , and snapshots from ab initio molecular dynamics (AIMD) simulations. A test set of structures was further generated to validate the generalizability of the fitted model, which is about 10% of that for the training set.
A comparison of the DFT-and MPEA SNAP-predicted energies and forces for both training and test sets is shown in Fig. 2. The corresponding MAE values in predicted energies and forces from the SNAP model relative to DFT for the elemental, binary, ternary, and quaternary systems are displayed in Supplementary Table 4. An excellent fit was obtained for the MPEA SNAP model, with a unity slope in both energies and forces with respect to DFT. The overall training and test MAE for energies and forces are within 6 meV/atom and 0.15 eV/Å, respectively. More critically, this high performance is achieved consistently across all sub-chemical systems, i.e., there is no obvious bias in performance for any particular chemistry.
The MPEA SNAP was further validated by computing various properties of the elements and multi-component systems, as presented in Table 1. While the MPEA SNAP model systematically overestimates the melting points for all elements, they are still in qualitative agreement with the experimental values. The elastic moduli predicted by the MPEA SNAP model are in extremely good agreement with the DFT for all four elemental systems, with errors of <10% except for the shear modulus of Mo (−11.3%). The MPEA SNAP also performs very well on the NbMoTaW SQS, except for a slight overestimate of c 44 and the corresponding shear modulus. It should be noted that only strained elemental structures, and not strained SQS structures, were included in the training data. Therefore, this excellent performance on the quaternary NbMo-TaW SQS is an important validation test for the generalizability of the MPEA SNAP model.  24,25 . α are the parameters (hyperparameters) needed for bispectrum calculations, e.g., the atomic weight S k atom for each element. β are the linear SNAP model parameters. ω are data weights from different data groups.  Generalized stacking fault (GSF) energies Metastable stacking faults play a critical role in the dissociation of dislocations in bcc metals 30,31 . The γ surfaces represent energies of GSFs, formed by shifting two halves of a crystal relative to each other along a crystallographic plane 32 . The MPEA SNAP model was used to compute a section along the {110}γ surface in the <111> direction (see Fig. 3a) for all four elemental bcc metals and the NbMoTaW SQS. Figure 3b shows the comparison between the MPEA SNAP results and previous DFT studies [33][34][35][36] . We can see that the overall agreement between SNAP and DFT results for all four elemental systems is excellent, with only a slight overestimation for the W. No local minima that would indicate the existence of metastable stacking faults was found in the elemental metals and the NbMoTaW SQS. The stacking fault energies of the NbMoTaW SQS are smaller than that of W and Mo and much larger than those of Ta and Nb.
Dislocation core structure One of the important characteristics of the screw dislocation in bcc metals is the core structure. Two types of 1/2 <111> screw dislocation core structures have been reported in previous calculations for bcc metals [36][37][38][39][40] , the degenerate core and the non-degenerate (or symmetric) core. A major discrepancy between ab initio methods and classical force fields is that the former predicts the non-degenerate configuration to be the core structure of the 1/2 <111> screw dislocation in bcc metals 36,39,40 , while the latter generally finds the degenerate core 37,38 . The MPEA SNAP model accurately predicts the non-degenerate core structure for the 1/2 <111> screw dislocation for all the bcc elements, consistent with DFT. Figure 4 shows classical differential displacement 41 plots at 0b and 4b along the dislocation line (b is the length of the Burgers vector) for the 1/2 <111> screw dislocation of the NbMoTaW SQS MPEA. It may be observed that  there is substantial variation in the core structure in the MPEA due to the different local chemical environments. The core structure is compact at 0b but shows more non-compact characteristics at 4b. Similar variations of the core structures in MPEA due to compositional fluctuations have been reported recently 42,43 . Nevertheless, sampling hundreds of local environments in the SQS cell indicate that the majority of local environments exhibit a compact core as shown in Fig. 4a, and only a very few local environments exhibit more non-compact characteristics as illustrated in Fig. 4b. These findings obtained from the SNAP calculations are consistent with the recent DFT studies in the same quaternary system 43 .
Critical resolved shear stress (CRSS) of screw and edge dislocations CRSS to move dislocations is closely related to the strength of the materials.  Table 2.
Generally, the MPEA SNAP model predicts a very high CRSS for screw dislocation in NbMoTaW SQS, much larger than that of Nb, Ta, and Mo and comparable with that of W. The most interesting observation, however, is that the MPEA SNAP CRSS for the edge dislocation in the NbMoTaW SQS is also much higher than those of the elemental components and is about~20% of the CRSS of the screw dislocation. This greatly reduced anisotropy between screw and edge mobility suggests that the edge dislocation may play a more important role in the deformation of the bcc MPEA compared to in bcc elements.

Segregation and SRO
The validated MPEA SNAP model was applied to long-time, largescale simulations of both single-crystal and polycrystalline models of the NbMoTaW MPEA. The single-crystal and polycrystal models were constructed using supercells of dimensions 15.5 × 15.5 × 15.5 nm (48 × 48 × 48 conventional cell) and 11 × 11 × 11 nm (Fig.  5a), respectively, with a randomized elemental distribution with equiatomic quantities, i.e., 25% each of Nb, Mo, W, and Ta (see Fig.  5b). Hybrid MC/MD simulations were then performed to obtain low-energy microstructures for the quaternary NbMoTaW MPEA at room temperature (see "Methods" for details). One important property that can be analyzed is the structural characteristics, such as pair correlation functions. For the singlecrystal MPEA, the partial pair correlation functions are plotted in Fig. 6a for both the random structure and the structure after equilibration in the MC/MD simulations. The dominant nearestneighbor correlations in the structure after MC/MD equilibration are between elements in different groups in the periodic table, with Ta-Mo being the highest, followed by Ta-W, Nb-Mo, and Nb-W; the correlations between elements within the same group (Ta-Nb and Mo-W) are much lower. This is consistent with the enthalpies of pairwise interactions (see Supplementary Table 5). The non-uniform correlations indicate the existence of local chemical order in the structure at room temperature. The computed pairwise multi-component SRO parameters 48,49 (see "Methods") are presented in Table 3. It was found that three intergroup elemental pairs-Ta-Mo, Ta-W, Nb-Mo-have large   Fig. 5c shows a snapshot of the polycrystalline model after equilibration in the MC/MD simulations. Clear segregation of Nb (blue atoms) to the GBs can be observed, while there is evidence of enrichment of W in bulk. Table 4 provides the atomic percentages for each element in the GBs and bulk before and after the MC/MD simulations. Starting with an initial equal distribution of approximately 25% for all elements in both GBs and bulk, the percentage of Nb in the GBs increases to~57.7%, while the percentage for W and Ta decreases to~2% and~13.6%, respectively. Correspondingly, the percentage for W and Ta in the bulk regions increases to~32% and~28%, respectively, while that for Nb decreases to~16%. The corresponding partial radial distribution functions are also plotted for this polycrystalline model after equilibration in the MC/MD simulations, as shown in Fig. 6b. We can clearly see a large decrease for the nearest-neighbor peak of Nb-W (blue curve) compared to the single-crystal model due to the segregation effects. In addition, the peaks for the Nb-Nb pair are broadened owing to the more disorder characteristics of Nb-Nb pairs segregated into the GBs. The segregation of W into the grain interior also leads to a large increase for the nearest-neighbor peak of W-W. The calculated SRO parameters (see Table 3) also indicate very different interactions between the polycrystal and single-crystal MPEAs. For example, α 1 NbÀW changes from a small negative value (−0.05) to a large positive value (0.34), due to the tendency of Nb and W to segregation into the GB and bulk regions, respectively.
MD simulations using the MPEA SNAP model were performed to generate uniaxial compressive stress-strain responses of nanocrystalline models of the elements as well as the random NbMoTaW MPEA and the equilibrated NbMoTaW MPEA after MC/ MD simulations, as shown in Fig. 7a. Among the elements, W has the highest strength, followed by Mo, and Ta and Nb being much weaker. The random MPEA has a strength that is substantially higher than that of Mo, Nb, and Ta. Most interestingly, the MC/MD-equilibrated NbMoTaW MPEA exhibits substantially higher strength than the random solid solution MPEA and close to that of W, the strongest elemental component. These results are consistent with previous experimental measurements 51,52 , which found that nanopillars of the MPEA has comparable compressive stress-strain curves with those of W at similar diameters (Fig. 7b).

DISCUSSION
We have developed a highly accurate SNAP for the four-component Nb-Mo-Ta-W system and applied it in large-scale, long-time simulations of both single-crystal and polycrystal NbMoTaW MPEAs. The accuracy of the MPEA SNAP model has been thoroughly evaluated based on not just accuracy in energy and force predictions but also in the reproduction of key mechanical properties, such as the elastic constants, dislocation core structure, and CRSS.
In the single-crystal MPEA, we find strong evidence of a reduced screw/edge anisotropy in the calculated CRSS. This finding supports recent observations that edge dislocations may play a far more important role in MPEAs as compared to bcc elements 53,54 . For example, Mompiou et al. 59 have found that edge dislocations are sluggish at room temperature in the Ti 50 Zr 25 Nb 25 alloy, indicating the comparable role in the strength of edge to screw dislocations. Maresca and Curtin 60 have also proposed that the random field of solutes in the highconcentration alloys has been found to create large energy barriers for thermally activated edge glide and established a theory of strengthening of edge dislocations in BCC alloys.
Large-scale, long-time simulations using the MPEA SNAP model have also provided critical new insights into the interplay between segregation, SRO, and mechanical properties of the polycrystalline NbMoTaW MPEA for the first time. First, it was found that there is a clear tendency for Nb to segregate to the GB, accompanied by enrichment of W in the bulk. Similar elemental segregation to GBs have also been observed in FeMnNiCoCr MPEA after aging heat treatment in a recent experiment 55 . This effect can be explained by considering the relative GB energies of the different elements. The current authors have previously developed a large public database of GB energies for the elemental metals using DFT computations 56 . As shown in Supplementary Fig. 3, Nb has the lowest GB energy and W has the highest among the four component elements. Hence, Nb segregation to the GB region and W enrichment in the bulk is driven by a thermodynamic driving force to lower the GB energies.
In turn, Nb segregation has a substantial effect on the observed SRO in the NbMoTaW MPEA. As can be seen from Table 3, the Nb-W SRO parameter changes from a small attractive interaction in the single crystal to a strong repulsive interaction in the polycrystal due to Nb segregation to the GB and W to the bulk. The SRO parameters of other pairs of elements are also intensified in magnitude. An increase in SRO has been found to lead to increased barriers to dislocation motion, leading to greater strength in the fcc NiCoCr MPEA 49 . Indeed, a similar effect is observed in the polycrystalline MPEA, where the equilibrated NbMoTaW MPEA with SRO exhibiting substantially higher strength than the random solid solution NbMoTaW MPEA, with a strength approaching that of W (Fig. 7a). The von Mises strain distribution 57,58 at a low applied strain of 3.0% is plotted in Fig. 8. It can be observed that the von Mises strain distribution is localized in the GB region for both the random solid solution and the equilibrated MPEA. However, the MC/MD-equilibrated polycrystal with Nb-rich GBs shows much smaller von Mises strains than the random solid solution. Similar GB stability-induced strengthening has also been observed experimentally in Ni-Mo nanograined crystals 59 .
To conclude, this work has highlighted that accurate treatment of MPEA chemistry at the inter-and intra-granular is necessary to reveal the subtle interactions between segregation and SRO and their consequent effect on mechanical properties. Interestingly, Nb enrichment of the GB coupled with and intensification of the SRO are predicted to the enhancement of the strength of the NbMoTaW MPEA over a random solid solution. These findings point to the potential for leveraging on both composition and processing modification to tune the GB composition and bulk SRO to tailor the mechanical properties of MPEAs.

Bispectrum and SNAP formalism
The SNAP model expresses the energies and forces of a collection of atoms as a function of the coefficients of the bispectrum of the atomic neighbor density function 60 . The atomic neighbor density function is given by: where δ(r − r ik ) is the Dirac delta function centered at each neighboring site k, the cutoff function f c ensures a smooth decay for the neighbor atomic density to zero at the cutoff radius R c , and the atomic weights S k atom distinguishes between different atom types. This atomic density function can be expanded as a generalized Fourier series in the four-dimensional hyperspherical harmonics U j m;m 0 as follows: where u j m;m 0 are the coefficients. The bispectrum coefficients are then given by: where C jm 0 j 1 m 0 1 j 2 m 0 2 are the Clebsch-Gordon coupling coefficients. In the linear SNAP formalism 22,24 , the energy E and force on atom j f j are expressed as a linear function of the K projected bispectrum components B k and their derivatives, as follows: where α is the chemical identity of atoms, N α is the total number of α atoms in the system, and β α,k are the coefficients in the linear SNAP model for type α atoms.
The key hyperparameters influencing model performance are the cutoff radius R c for bispectrum computation, atomic weight S k atom for element k (Nb, Mo, Ta, or W), and the order of the bispectrum coefficients j max . In this work, the LAMMPS package 61 was used to calculate the bispectrum coefficients (the features) for all the training structures 22 . An order of three for the bispectrum coefficients (j max = 3) was used, consistent with previous works 22,24,25,60 . The cutoff radius R c and atomic weight S k atom were optimized during the training of the model.

Training data generation
One critical factor for developing an effective and robust potential is a diverse training data encompassing a good range of atomic local environments. For a quaternary potential, the training data should include the elemental, binary, ternary, and quaternary systems. The detailed structure generation for each system is provided as follows. Fig. 6 Partial radial distribution function. The plots of pair correlation functions g αβ at T = 300 K for a a single-crystal NbMoTaW MPEA and b a polycrystalline NbMoTaW MPEA. Table 3. Pairwise chemical short-range order parameters α 1 (see "Methods") for the NbMoTaW MPEA single crystal and polycrystal after MC/MD equilibration at room temperature.

Pairs
Single crystal Polycrystal  • Snapshots from NVT AIMD simulations of the bulk 3 × 3 × 3 supercell at room temperature, medium temperature (below melting point), and high temperature (above melting point). In addition, snapshots were also obtained from NVT AIMD simulations at room temperature at 90% and 110% of the equilibrium 0 K volume. Forty snapshots were extracted from each AIMD simulation at intervals of 0.1 ps; 2. Binary systems (Nb-Mo, Nb-Ta, Nb-W, Mo-Ta, Mo-W, Ta-W) • Solid solution structures constructed by partial substitution of 2 × 2 × 2 bulk supercells of one element with the other element.
Compositions of the form A x B 1−x were generated with x ranging from 0 to 100 at% at intervals of 6.25 at%.
For the binary solid solution structures with each doping percentage, we performed a structure relaxation for all symmetrically distinct structures. Both the unrelaxed and relaxed structures were included in our data set. For the ternary and quaternary systems, our data set also includes structures generated by permuting the elements in the generated SQS, as well as the relaxed structures (including the intermediate structures) by optimizing the generated SQS.
The test set of structures was generated by extracting additional snapshots from all previous AIMD simulations at intervals of 0.1 ps. We also generated additional binary solid solution structures by partial substitution of one element with the other element in a 2 × 2 × 1 supercell. The substitution percentage ranges from 0 to 100 at% at intervals of 25 at%. The total number of test structures is about 10% of that for training data.

DFT calculations
We performed the DFT calculations using the Perdew-Burke-Ernzerhof 66 exchange-correlation functional and projector-augmented plane wave 67 potentials as implemented in the Vienne ab initio simulation package 68 . The plane-wave cutoff energy was 520 eV, and the k-point density was 4 × 4 × 4 for 3 × 3 × 3 supercells. The energy threshold for self-consistency and the force threshold for structure relaxation were 10 −5 eV and 0.02 eV/Å, respectively. A single Γ k point was applied for non-spin-polarized AIMD simulations. However, we used the same parameters as the rest of the data for the energy and force calculations on the snapshots. The Python Materials Genomics (pymatgen) 69 library was used for all structure manipulations and analysis of DFT computations. Fireworks software 70 was applied for the automation of calculations.

SNAP model fitting
The fitting workflow for the quaternary SNAP model is illustrated in Fig. 1, in which we adopt the potential fitting workflow for elemental SNAP model developed in 24 as an optimization unit. This optimization unit contains two optimization loops. The inner loop optimizes the ML model parameters (β in Fig. 1) by mapping the descriptors (bispectrum coefficients) to DFT-calculated formation energies and forces. The formation energies are defined as, E form = E TOT − ∑ el = Nb,Mo,Ta,W N el E el , where E TOT is DFT-calculated total energy of the system, N el is the number of atoms in the system for the each type of element, and E el is the energy per atom in the corresponding elemental bulk system. The hyperparameters are optimized in the outer loop by minimizing the difference between the model-predicted material properties, i.e., elastic tensors, and the DFT-computed values. These hyperparameters include the data weight (ω in Fig. 1) from different data groups and the parameters (α in Fig. 1) used in bispectrum calculations, i.e., the radius cutoff R c , and atomic weight S k atom . The fitting algorithm for each loop is the same as previous  works 24,25 with the least-squares algorithm for inner loop and the differential evolution algorithm 71 for the outer loop. For the quaternary NbMoTaW alloy system, there are eight hyperparameters (R Nb c , R Mo c , R Ta c , R W c , S Nb atom , S Mo atom , S Ta atom , S W atom ) in the bispectrum calculations, two for each element. A more efficient step-wise optimization was performed. In the first step, we performed a series of independent optimization of the radius cutoff R c for each elemental SNAP model, i.e., Nb, Mo, Ta, and W. The optimized radius cutoffs are then used as the radius cutoff for the quaternary NbMoTaW SNAP model. In the second step, a grid search was performed for the atomic weight for each element. We initially fix the data weight according to the number of data points for each data group and perform a quick grid search for the atomic weight of each element by only running the inner loop. The grid range is confined between 0 and 1.0 with an interval of 0.1 for each atomic weight. We then select the first ten combinations of the atomic weights of the four elements (see Supplementary Table 2) with the best accuracy in energy and force predictions to conduct a full optimization, including the optimization of the data weight in the outer loop. The best model (with the highest accuracy in energies and forces) was selected from the ten fully optimized models.

Atomistic simulations
Atomistic simulations using the MPEA SNAP model were performed using the LAMMPS code 61 . Specifically, • Melting points. The solid-liquid coexistence approach 72 was used for melting temperature calculations. We use a 30 × 15 × 15 bcc (13,500 atoms) supercell for each system to perform MD simulations under zero pressure at various temperatures. One fs was set for the time step. We performed the MD simulations for at least 100 ps to ensure the correct conclusion. The temperature, at which the initial solid and liquid phases were at equilibrium without interface motion, was identified as the melting point.
• GSF energies. GSF energies were performed using a large supercell containing about 36,000 atoms. The supercell was set to be periodic along <111> and <112> directions in the {110} plane and non-periodic along <110> direction.
• Dislocation core structure. To study dislocation core structure and dynamics, we inserted a 1/2 [111] screw dislocation with line direction z = [111], glide direction x = [11−2], and glide plane normal y = [−110] into a cylinder supercell with a radius 10 nm. The length of the dislocation line of one periodicity along <111> direction is 8b (b is the Burgers vector) in the supercell. The quaternary cylinder supercell is constructed from a 54-atom SQS. The dislocation was inserted by deforming the atomic positions according to the linear elasticity theory. Rigid boundary conditions were used by creating a layer of atoms fixed in their unrelaxed position outside of the inner cylinder region with radius 9 nm of mobile atoms. This method with this configuration has been used to study dislocation properties in previous works for bcc metals 34,73 . Similarly, a 1/2 [111] edge dislocation with the line direction along z = [11−2] could be introduced inside the cylinder supercell. Energy minimization was performed using periodic boundary conditions along the dislocation line direction (z direction) and fixed boundary conditions along the other two directions (x and y directions). To measure the CRSS or Peierls stress for dislocation motion at T = 0 K, we applied increasing homogeneous shear strain in small increments and determined the stress value at which the dislocation moves from its initial position as the CRSS. For MPEA, we record the largest shear stress within one periodicity.
• Polycrystal simulations. The initial polycrystal model was generated using the Voronoi tessellation method 74 implemented in the Atomsk code 75 . We constructed a 11 × 11 × 11 nm supercell and randomly inserted six GBs with an average grain diameter of about 7.5 nm. Periodic boundary conditions were imposed in all directions. At the GBs, pairs of atoms with a distance <1.5 Å were removed. A hybrid MC/MD simulation was performed on the MPEA polycrystal for 1.5 ns. The MD simulations were carried out in NPT ensembles at 300 K and zero pressure. The time step was set to 5 fs. In the MC run, the sample was initially generated by randomly inserting the four elements into the supercell. Trial swaps were performed by randomly selecting an atom and replacing it with another species. The relative mole fractions of different elements are kept constant. The trial moves are accepted or rejected based on the Metropolis algorithm. A swap attempt was conducted at every MD step.
• Stress-strain simulations. Uniaxial compressive deformation simulations were carried out with a strain rate 5 × 10 8 s −1 . The time step was set to 1 fs, and simulations were performed under NPT ensemble at room temperature.

Chemical short-range-order parameters
The definition of the pairwise multicomponent SRO parameter is α k ij ¼ ðp k ij À c j Þ=ðδ ij À c j Þ; (6) where k denotes the kth nearest-neighbor shell of the central atom i, p k ij is the average probability of finding a j-type atom around an i-type atom in the kth shell, c j is the average concentration of j-type atom in the system, and δ ij is the Kronecker delta function. For pairs of the same species (i.e., i = j), a positive α k ij means the tendency of attraction in the kth shell and a negative value suggests the tendency of repulsion. For pairs of different elements (i.e., i ≠ j), it is the opposite. A negative α k ij suggests the tendency of j-type clustering in the kth shell of an i-type atom, while a positive value means the repulsion.

DATA AVAILABILITY
To ensure the reproducibility and use of the models developed in this paper, all data (structures, energies, forces, etc.) used in model development as well as the final fitted model coefficients have been published in an open repository (https://github. com/materialsvirtuallab/snap).

CODE AVAILABILITY
The DFT calculations were performed with the Vienna ab initio simulation package. The LAMMPS package was used to perform MD/MC simulations. All the other codes that support the findings of this study are available from X.-G.L. (xil110@ucsd.edu) upon reasonable request.