Ab initio framework for deciphering trade-off relationships in multi-component alloys

While first-principles methods have been successfully applied to characterize individual properties of multi-principal element alloys (MPEA), their use in searching for optimal trade-offs between competing properties is hampered by high computational demands. In this work, we present a framework to explore Pareto-optimal compositions by integrating advanced ab initio-based techniques into a Bayesian multi-objective optimization workflow, complemented by a simple analytical model providing straightforward analysis of trends. We benchmark the framework by applying it to solid solution strengthening and ductility of refractory MPEAs, with the parameters of the strengthening and ductility models being efficiently computed using a combination of the coherent-potential approximation method, accounting for finite-temperature effects, and actively-learned moment-tensor potentials parameterized with ab initio data. Properties obtained from ab initio calculations are subsequently used to extend predictions of all relevant material properties to a large class of refractory alloys with the help of the analytical model validated by the data and relying on a few element-specific parameters and universal functions that describe bonding between elements. Our findings offer crucial insights into the traditional strength-vs-ductility dilemma of refractory MPEAs. The proposed framework is versatile and can be extended to other materials and properties of interest, enabling a predictive and tractable high-throughput screening of Pareto-optimal MPEAs over the entire composition space.

However, certain alloy compositions have been identified to have a more favorable combination of strength and ductility than others [12].For instance, the ductilisation effect of rhenium (Re) and iridium (Ir) in tungsten is a well-known example that keeps the overall strength of the material virtually unchanged [13][14][15].Consequently, the quest for discovering new compositions of refractory MPEAs with improved mechanical properties remains an important challenge in materials science.The conventional trial-and-error method of materials design is too laborious for complex MPEAs.This is where computerguided design opens new opportunities for efficient exploration of the vast composition space of these materials.
In recent years, many works have introduced computational methods to identify and characterize alloys with tailored properties [16][17][18][19][20][21][22][23], such as the solid solution strengthening, hardness, ductility index, and thermodynamic stability.Moreover, some of these methods have been employed in the context of the optimization of several target properties [16,[24][25][26], with the relevant material properties being obtained from various modeling approaches using experimental and theoretical input.For example, methods based on machine learning (ML) use various sets of descriptors partly derived from ab initio calculations in order to fit experimental data [27][28][29].This approach, however, often requires large data sets and, most crucially, it is intrinsically interpolative and cannot give predictions outside the domain covered by an existing data set.A more transferable methodology is to use properties of elements, or simple compounds, for extrapolation of target properties to multi-component properties, mainly using the rule-of-mixture [30][31][32].Despite its appeal due to its simplicity, this approach cannot capture the full complexity of a system in the case of intricate interactions (e.g., due to abrupt changes in the near-Fermi level band structure [33]).
A typical method of choice for truly predictive calculations of alloy properties is density functional theory (DFT).However, explicit evaluation of fracture or yield strength in alloys is mostly out of reach for available DFT implementations.Although this may become possible with general-purpose interatomic potentials [34,35], a fairly robust description can be obtained by using appropriate models, which have been shown to accurately predict strengthening [36][37][38] and ductility [39][40][41][42].These models combine linear elasticity and microscopic principles, and the only required input parameters are cer-tain material properties.Specifically, apart from basic equilibrium properties, such as the molar volume and elastic constants, the strengthening model requires misfit volumes related to concentration derivatives of the alloy volume, while the ductility model relies on surface and unstable stacking fault (USF) energies computed for specific orientations.Importantly, all of these parameters are accessible through simulations using first-principles methods, yet, such simulations are still challenging for disordered MPEAs.
The most common approach to ab initio alloy modeling is to represent a disordered alloy using special quasirandom structures (SQS) or similar supercell representations that try to mimic perfect randomness [43][44][45].However, this method has several inherent limitations for MPEAs.First, the computational effort grows exponentially as the number of alloy components increases.Moreover, properties influenced by global symmetry or local chemical environment [46,47], such as surface energies or elastic constants, present significant challenges for convergence and necessitate the use of even larger cells, limiting the number of components in the best case to three or four.Calculating properties for arbitrary compositions, especially when components are present in low concentrations, is even more cumbersome.
One alternative methodology is the coherent potential approximation (CPA) [48,49], whose advantages are its extremely small computational expenses, irrespective of the number of components, as well as its ability to treat arbitrary alloy concentrations, providing a straightforward way for accurate evaluation of concentration derivatives.Unfortunately, the single-site nature of CPA renders it impossible to take into account effects of local atomic relaxations that are crucial for a correct characterization of defects, such as USFs, or inhomogeneous structures, e.g., surfaces.
Another alternative to traditional supercell methods are machine-learning interatomic potentials (MLIPs) that are trained on ab initio data [50][51][52][53].The size of supercells is not limited with MLIPs, which is especially important when simulating extended defects like interfaces, where proper configurational averaging becomes an issue with MPEAs.In contrast to the previously proposed way of passively training general-purpose interatomic potentials on a pre-defined fixed training set (e.g., [34,35,54,55]), whose accuracy is not sufficient for our goals, especially for parts of the configurational space not covered by the training set, the use of active learning enables automatic generation of MLIP training sets specifically tailored to a given class of systems and a given set of properties.This allows one to predict any material property of interest with DFT accuracy over the entire composition space of an alloy [39,56] with a minimal amount of DFT calculations.
In this work, we combine the two approaches-CPA and actively learned Moment Tensor Potentials (MTP), a class of MLIPs-to characterize solid solution strengthening and ductility of refractory MPEAs.In particular, we use CPA to compute the elastic constants and misfit volumes, needed for computing the critical resolved shear stress (CRSS) within the model of Maresca and Curtin [36,37].MTPs are used to simulate large supercells, necessary to calculate surface and USF energies which are needed to parameterize a model for intrinsic ductility based on the Rice-Thomson theory [40].
The crucial point of our approach is to enhance the simulation efficiency by employing the most suitable method for different sets of properties.Specifically, CPA's quantum mechanical treatment of solid solutions and its analytical nature eliminates finite-size cell effects and sampling errors related to random atomic configurations, making it excel in calculations of the concentration dependence of properties that are relatively insensitive to local relaxations, such as elastic constants and misfit volumes.In comparison to supercell calculations with MLIPs, CPA is especially superior in accuracy when dealing with abrupt changes of the Fermi-surface topology, which is often encountered in bcc alloys and can have a significant impact on the concentration dependence of properties.On the other hand, ab initio based MTP methodology surpasses CPA in dealing with defects and special geometries, such as surfaces, due to its capability to incorporate local atomic relaxations.
Another important part of our framework is a simple model that is only slightly more complicated than the rule of mixtures but can accurately capture the concentration dependence of all the quantities relevant for strengthening and ductility over the entire composition space of a large class of refractory alloys.The model is based on simple element-specific descriptors and universal functions of the d-electron valence, with a few adjustable parameters that can be easily fitted to data from simple compounds.Importantly, the model is based on fundamental properties of bonding in transition metals, which makes it capable of extrapolative predictions, rendering it as a potentially useful tool for exploring the entire space of refractory MPEAs.Moreover, the analysis of higher-order contributions beyond the first order can provide an insight into interactions between individual alloy components.
As a proof of concept for the flexibility of such an approach, we perform multi-objective optimization (MOO) of strength and ductility for three-and four-component alloys.Then, we use these results to validate the model and extend our predictions to a larger class of alloys.By analyzing the results, we reveal several potential pitfalls of conventional methods and show how they are easily overcome within our approach.

A. Multi-objective optimization of strength and ductility
We consider homogeneous alloys whose strength can be characterized by the CRSS due to impeded mobil-ity of dislocations interacting with the alloy components that act as a random field of solutes.Within the Maresca-Curtin model [37], the CRSS is given by τ y ≡ τ y (V, {∆V n }, C ij ), where V is the specific volume, ∆V n are the misfit volumes for each component n, and C ij are the elastic constants (see Methods IV B for details).If the material parameters are obtained from first principles, the model allows for parameter-free predictions of yield strength without relying on experimental input.
The strength of a solid solution is known to compete with its ductility.The ductility of a refractory MPEA is mainly determined by the intrinsic fracture behavior of single-crystalline grains.The ability of a material to undergo plastic deformation affects the way unavoidable defects grow to cracks and propagate through the material.If a crack is blunting (propagating) upon external loading, the material is said to be intrinsically brittle.If, instead, a dislocation is emitted from the crack tip, the material is said to be intrinsically ductile.This competition between blunting and dislocation emission on active slip systems in bcc alloys can be quantified using the Rice-Thomson model [40,57].Within the Rice-Thomson model, ductility is described by means of a ductility index, D, which depends on the ratio of the unstable stacking fault energy γ USF and the surface energy γ surf , times a prefactor depending on the crack orientation and the elastic constants.In this theory, lower D values indicate increased alloy ductility, and D < 1 signifies the absence of intrinsic brittleness.A recent application of the ductility index by Mak et al. [41] has led to the development of an RT ductility criterion for refractory metal alloys.Their study showcased the correlation of D with measured fracture strain and its predictive capability for ductility trends.
Improved ductility is expected for alloys with lower stacking fault energies, which facilitates dislocation emission.At the same time, the opposite is true for improved strength, and to investigate this competition, we leverage a Bayesian multi-objective optimization framework to uncover sets of Pareto-optimal compositions, with the objective functions obtained from first principles.The flowchart in Fig. 1 presents an overview of the workflow developed in this work.The two main ingredients of the workflow are 1) the calculations of the objectives, the CRSS, τ y , and the ductility index, D (Figure 1a), and 2) a Bayesian optimization loop based on the nondominated sorting genetic algorithm (NSGA-II) [58] for finding the Pareto-set approximation (Figure 1b,c); see Methods IV D for details.The key feature here is that NSGA-II utilizes an intermediate surrogate model (Gaussian process in our case), significantly reducing the need for direct evaluations of the objective functions.
We use CPA to calculate the elastic constants and volumes because these quantities are insensitive to local atomic relaxations, and evaluation of concentration derivatives is much more efficient with this approach compared to SQS calculations.Furthermore, a similar methodology has already been successfully applied to computing CRSS of iron-group MPEAs [59,60].
On the other hand, the planar defects are strongly affected by local lattice relaxations, and we, therefore, resort to a supercell approach for these quantities.However, carrying out high-throughput calculations with the conventional SQS is computationally unwieldy, especially for properties requiring non-trivial structure relaxations (e.g., surfaces).Therefore, we use actively-learned MTPs that are able to predict material properties with near-DFT accuracy [39,56], requiring only a small fraction of the computational expenses of a direct DFT calculation.This allows us to perform structure relaxation on supercells with hundreds of thousands of atoms at a cost that is only slightly higher than that of classical interatomic potentials, such as the embedded atom method.In comparison to previously published works [27,41,61,62], our approach consistently produces comparable USF energies for specific alloy compositions (see Results and Discussion) but we can also treat any other composition of the alloy once the MLIPs are trained.Moreover, we anticipate greater accuracy for systems affected by Fermi surface nesting ( [63]) or for alloys experiencing local lattice reconstructions ( [64]) because our methodology does not impose constraints on cell size or lattice relaxations.Indeed, in case of surface energies, the differences between MTP-based and SQS calculations become much more pronounced (sometimes, more than 200 mJ/m 2 , which is more than 10% of the typical value of γ surf ), mainly due to a considerably slower convergence of the surface energy with the cell size, as will be discussed in more details later on (see Section III A).

B. Search for Pareto-optimal alloys
We now combine the calculations of the CRSS and the ductility index, D, and perform a search for Paretooptimal alloys in the composition space, according to the workflow shown in Fig. 1.We first analyze the MoNbTi alloy, whose search (design) space can be easily visualized using ternary diagrams.We constrain the search space to a concentration range of 0.0 ≤ c i ≤ 0.75 for each component i.This is motivated by the fact that pure Ti is stable in the hcp structure (α-Ti) under ambient conditions and it has limited miscibility with the other two elements, Mo [67] and Nb [68] (for more detailed analysis of phase stability, see Supplemental Materials -SM S1).Moreover, dilute alloys offer minimal potential for solid solution strengthening.
The constraints are handled by penalty functions within NSGA-II, which penalize undesirable solutions by reducing their fitness values in proportion to the degree of constraint violation.In principle, this approach could be extended to define more complex regions in the search space that should be avoided, e.g., regions of thermody- namic instability or other unfavourable properties.
Fig. 2a shows the sampling points in the search space of MoNbTi, alongside with the Pareto set approximation (PSA), which is a set of compositions that are nondominated by each other but are superior to the other compositions in the Pareto sense.This can be seen more clearly in Fig. 2b, where we show the Pareto front approximation (PFA) at the bottom of the feasible region, resulting from the sampling in the objective space.There is a continuous mapping between PFA and PSA.The PFA's left portion corresponds to PSA's right portion, and vice versa.On the right side of PSA and the left side of PFA are stronger but less ductile alloys.In both PFA and PSA, the color-coded scheme designates yellow for high strength and pink for high ductility.Generally, we observe that relative changes in the ductility index across the concentration space are not as severe as changes in the strength.The lower limit of the strength for Paretooptimal alloys within the constraints given above is about 50 MPa, while the strongest alloys have τ y ≈ 260 MPa.On the other hand, the corresponding ductility index, D, changes from 1.23 to 1.35 (with the error of D estimated to be less than 0.015, see SM S3 E).It is also worth noting that the Pareto set avoids the region around the equimolar composition.
Regions of higher strength can be found close to the equimolar MoNb alloy with small additions of Ti.Reducing the Mo content in the alloy enhances its ductility and reduces the likelihood of brittle fracture.On the other hand, the alloy must contain a certain amount of Mo to ensure a reasonable strengthening effect.More ductile alloys can be found in a region close to the Nb 0.6 Ti 0.4 binary.The improved ductility comes at a price of the strength that is almost five times smaller than the maximum value.
Pareto-optimal alloy compositions in the threecomponent MoNbTi alloy system exhibit a continuous and smooth transition across the entire design space, forming a connected Pareto set.To investigate nontrivial alloying effects, we expand the composition space by introducing Ta, anticipating Fermi surface nesting effects caused by Ta alloying.Similarly to the case of MoNbTi, we confine the search space to concentration values ranging from 0.0 to 0.75 for each component, with the exception of Ti, which is restricted to a maximum atomic fraction of 0.55.This limitation is motivated by the limited solubility of Ti, as described in Ref. [69] (see also SM S1).The results of the optimization, represented by the PFA and sampling points, as illustrated in Fig. 3a reveal a more intricate concentration dependence of the objectives, which can also be seen in the t-SNE projections displayed in Figs.3b,c (see Methods IV E for details).
These plots indicate that, while the ductility index, D, is, again, minimal near binary NbTi, the strength, τ y , is maximal in two regions: 1) close to binary MoNb, with small additions of Ti, as in the MoNbTi case; 2) the vicinity of binary MoTa, which exhibits considerably higher strength than the 3-component alloy.This leads to two effectively disjoint Pareto subsets in the composition space of the MoNbTiTa alloy.
Indeed, clustering analysis has detected that parts of the Pareto front denoted by PFA1 and PFA2 in Fig. 3a can be associated with two separate regions in the concentration space.Despite a seemingly continuous crossover between these two regions in the objective space, the regions are truly disjoint in the composition space.The right, ductility-favoring portion of the Pareto front, PFA1, corresponds to alloys similar to the threecomponent MoNbTi system with minor additions of Ta.The optimal compositions in this region are characterized by the Nb content ranging from 0.4 to 0.6, by the Ti content ranging from 0.2 to 0.4, and moderate Mo concentration maintaining the strengthening.At the same time, the left portion of the Pareto front, PFA2, represents a separate region comprising alloys with high Mo (0.4−0.6) and moderate Ta (0.2−0.35) content, whereby Ta is responsible for a strength above 250 MPa, i.e., beyond those reachable by the MoNbTi system.An important feature of the PFA2 is a small knee corresponding to τ y ≈ 270 MPa.Starting from this point, one can either improve ductility without reducing the strength much, or enhance the strength without compromising on ductility.
In Figs.2a and 3a, we highlight three alloys from Ref. [65,66] which were mechanically tested.Both the equimolar MoNbTiTa and MoNbTi exhibit similar strength in experimental measurements, with the fourcomponent alloy showing slightly higher strength.The Nb-rich Mo 0.12 Nb 0.7 Ti 0.18 alloy demonstrates approximately 2/3 of the strength of the others.This behavior is also observed in our results.Additionally, ductility was indirectly assessed in these experimental works by measuring strain at fracture.The Nb-rich alloy exhibited the highest ductility, while the other two alloys showed lower ductility, yet remained non-brittle.Our ductility index corroborates this finding.All three alloys are lo-cated near the Pareto-front, with the Nb-rich alloy being located precisely on the ductile portion of the Pareto front.For MoNbTiTa, our results suggest that increasing the Mo and Nb content, while adjusting the Ta content, can improve the ductility of the equimolar MoNbTiTa alloy, without sacrificing the strength.

C. Virtual bond approximation model
So far, the main idea of the workflow was to get accurate and predictive results with high computational efficiency.Still, characterizing multiple alloy systems in this ways would require a lot of DFT calculations (MTP training and CPA).Instead, we would like to use the results obtained thus far for MoNbTiTa to analyze the strengthening and ductility behavior in a broader range of refractory alloys.One way to achieve this would be to define a set of "good" descriptors and train a black box ML model, which can then be used to obtain properties at arbitrary compositions.Such descriptors may encompass purely local structural characteristics [75][76][77], effective interatomic interactions [78] or local electronic structurerelated properties like moments of the local density of states and energies of simpler substructures [27,79].One significant disadvantage of such data-driven models is that they often lack interpretability and they are also only interpolative, meaning that one has to have a large amount of data covering enough of alloy systems to be  able to predict properties of an alloy not included in the initial data set.Physical models, on the other hand, are more transferable because they require fewer data points for training and provide clear insights into system behavior, fostering a deeper understanding of cause-and-effect relationships.Moreover, such models based on fundamental principles can be also considered as extrapolative, which prevents implausible outcomes, ensuring reliability.(a) Comparison of (a) C44 and (b) ∆V misfit volumes obtained from the rule-of-mixture (R), virtual bond approximation (M), and direct calculations (C).
One of the simplest examples of a physics-motivated model is a correlation between the valence electron count (VEC) and the ductility.In particular, a lower VEC is known to be associated with an enhanced intrinsic ductility due to a greater elastic softness [11,41,42,80].Brittle alloys are prevalent at a VEC of around 5 to 6, with an improved ductility beyond this range [80].A similar correlation can be seen for the shear modulus, G, as shown in Fig. 4b.The shear modulus increases more significantly than the bulk modulus, resulting in a larger shear-to-bulk modulus ratio (G/B) at higher values of VEC.This trend is frequently associated with embrittlement.
Despite its clear physical meaning as the bonding state filling factor, the VEC does not take into account interactions between alloy components.This, in particular, can explain why it fails to reveal any particular trend in misfit volumes (Fig. 4a) and it is, thus, not a useful descriptor for the CRSS, as can be seen in Fig. 4c.Another popular model -the rule of mixtures (ROM), also known as Vegard's law when applied to the lattice constant or volume -shares the same drawback (see, e.g., Fig. 5).
Our approach, novel in the context of strengthening and ductility, is based on a slightly modified version of the virtual bond approximation (VBA) [81] that can be considered as a generalization of classical bonding theory of transition metals [82] to alloys.The key idea of the VBA is to assign universal functions to local bonds between alloy components.These so-called virtual bond functions depend on the average d-valence of elements forming the bond, as well as element-specific bandwidths (equivalently, second moments of the density of states) of the corresponding elemental compounds.To second order in inter-component interactions (possible extension to higher orders is straightforward), the VBA takes the form where ξ({c i }) stands for B, C ′ , C 44 , V , or ∆E (for USF and surface energies), with c i being atomic fractions of components and ξ (k) representing the corresponding virtual bond functions of the average d-valences, N (k) ≡ {N i , N ij }, and of bandwidth parameters w (k) ≡ {w i , w ij } (see Methods IV F for details).Note that N i and w i are element-specific parameters, while the second-order parameters are defined as The unknowns ξ (k) are obtained using linear regression over a set of values of the associated quantity computed for a series of systems including pure elements (Cr, Mo, Nb, Ta, Ti, W, Zr) and equimolar binary alloys in the bcc structure.Since the inclusion of Ta causes deviations and abrupt changes in material properties attributed to Fermi surface transitions (e.g., see Fig. S4 in SM), we also calculate several three-and four-component alloys containing Ta. Significant deviations from both directly calculated and ROM values, serve as a reliable indicator of whether additional alloys should be included.This procedure was only deemed necessary for elastic properties.The fitted model is applicable to any alloy composed of the seven mentioned elements.
A comparable approach was presented by Ferrari et al. [20,83] where energies and material properties were fitted using polynomial functions of concentrations.A notable advantage of our method lies in the fact that we do not require parameters to be fit individually for each alloying system.Our approach promotes the coefficients from system-specific tabulated values to transferable parameters that encapsulate local bonding characteristics.For example, the VBA applied to iron-group MPEAs produced good results despite major complications due to the complex magnetic behavior of this class of alloys [59].
To demonstrate the performance of our VBA model, we compare direct and VBA calculations for the most important properties for six previously considered equimolar alloys: NbTiZr, NbTaTi, MoNbTi, CrMoTaTi, MoNbTaW and MoTaW.The results are presented in Fig. 5, where we compare calculated (labeled as "C"), VBA model ("M"), alongside ROM ("R") values for the misfit volumes, C 44 , and γ surf .The misfit volumes for equimolar compositions (Fig. 5b) are reproduced well by both the VBA model and the ROM, with the former giving systematically better values.However, the true advantage of using a model that includes inter-component interactions, becomes apparent when the composition is varied: In this case, the ROM often fails to reproduce the concentration trends, while the VBA model produces almost exact results corresponding to direct calculations (e.g., see Fig. S11 in SM).
Obtaining virtual bond functions, ξ (k) , for surface and USF energies is more difficult due to relaxation effects.To take them into account, we apply an additional rescaling of the functions, as described in Methods IV F. Test results can be found in SM S4 (Figs.S13 and S12).
To validate our VBA model against our ab initio results, we use the compositions previously generated during the optimization process for the MoNbTi and MoNbTiTa alloy systems, compute τ y and D using the material parameters obtained from the VBA model, and compare the results with direct CPA+MTP calculations.As can be seen from Fig. 6, the agreement is remarkably good for such a simple model, with minor deviations due to its inability to replicate strong non-linear effects in the properties upon changing the concentration.A more detailed analysis reveals that most of the discrepancies (especially for τ y in Fig. 6) are associated with Ta causing substantial non-linearity in the concentration dependence of the elastic properties (see SM S2 B).Such a non-linear behavior is due to Fermi-surface effects and pose a challenge to any simple model.
The validation on MoNbTiTa has shown the ability of the VBA model to to describe trends in multicomponent alloys based on the data from elemental and simple alloy compounds.However, a much stricter test (also addressing an eventual overfitting issue) would be to check if the model can predict properties for alloys containing elements not included in the training set.With this in mind, Hf and V (or any alloy containing them) have been excluded from the training set used to fit parameters ξ (k) .To evaluate properties for alloys involving these elements, the only input needed is the corresponding bandwidth parameter (corrected for the row in the periodic table) and the number of d-electrons, which are element-specific constants.
The model is then applied to a series of alloys containing nine d elements from three groups: Ti, V, and Cr.The results for the CRSS are shown in Fig. 7 where they are plotted against experimentally measured τ y [84][85][86][87] (tabulated data can be found in Table S2 of SM S4 B).One can see that the model captures trends very well in a large range of V-containing alloys.Some alloys display underestimated predictions, but the overall strength ranking remains accurate, reflected in a high Pearson coefficient.The overall underestimation is anticipated as the strengthening model gives the lower bound, as discussed in more details in SM S2 A (the effect can also be seen in Fig. 4).
The predicted ductility index is much more difficult to validate, as it is not a measurable quantity.Therefore, we compare it against the fracture strain for a range of alloys from Singh et al. [23].As one can see from the results in Fig. 8 (tabulated data can be found in Table S1 of SM S4 B), there is a significant anti-correlation between the measured strain at fracture and D, with the Pearson coefficient of -0.67.In critically assessing this outcome, one should keep in mind that it captures not only errors of the VBA model itself but also the qualitative nature of the ductility index, which is not directly related to the fracture strain, but rather, it is an indicator for the tendency to embrittlement.Furthermore, the D-criterion does not apply to intergranular fracture that can be facilitated by elements segregating at grain boundaries and reducing their cohesive strength.Nevertheless, it is important that most of the alloys on the brittle and ductile ends (having extreme values of D/fracture strain) are ranked correctly.Note also, that our results are very similar to an analogous plot of the fracture strain versus the ductility index in Singh et al. [23], if one takes into account that the ductility index in the reference is related to the inverse of our D. Overall, the prediction results suggest that the VBA is sufficiently accurate to be used for a qualitative assessment and analysis of ductility-strength trade-offs for refractory MPEAs, as will be discussed in the next section in more details.

A. Comparison to SQS results
Although our results on the CRSS and the ductility index are qualitatively consistent with previously published calculations, there are quantitative differences, which we now analyze by examining individual quantities involved in the models.
The key parameters for calculating the ductility index, D, are the energies of USF (for two orientations, {110} and {112}) and the surface ({100} and {110}).In this context, we compare the results of equimolar MoNbTi of our calculations with actively-learned MTPs to respective values from previous studies in Table I.
In particular, the USF are often used in theories predicting the mechanical properties of bcc metallic alloys, and its energy is therefore often studied.Our MTP approach yields values of 730 mJ/m 2 and 853 mJ/m 2 , for {110} and {112} respectively.The values for both orientations converge quickly with cell size, which validates the use of smaller cells in SQS calculations.
The SQS-based references values of Xu et al. [62] and values from MTPs trained on the large data set of Zheng et al. [61] render similar values, with rather small difference up to 40 mJ/m 2 .Errors on this scale can be attributed to the choice of a specific lattice parameter, statistical sampling errors and method-specific numerical errors (Supplemental Materials S3 D contain a detailed analysis of errors for surface and USF energies).Only Mak et al. [41] shows consistently higher values for all orientations.Notably, our analysis reveals the origin of the discrepancies among the reference values reported in the literature.Xu et al. [62] and Zheng et al. [61] calculated USFs by keeping the cell cubic with the lattice parameter of the alloy.In contrast, Mak et al. [41] relaxed the out-of-plane lattice vector to zero stress, along with the ionic positions perpendicular to the stacking fault, while keeping the fault plane lattice vector fixed.Replicating this procedure with our MTP on similar cell sizes revealed a systematic increase of approximately 70-100 mJ/m 2 in γ USF for various compositions.In comparison, alternative approaches, such as the surrogate model of Hu et al. [27], give comparable, but overestimated, values (the last row in Table I).Additionally, our calculations highlight the substantial contribution of local ionic relaxations to the USF energy, causing an energy change of about 25%.Despite this, neglecting local relaxations does not affect overall trends and relative changes upon concentration variations, potentially justifying the use of CPA calculations for qualitative studies (more detailed analysis provided in SM S3 A).In contrast to USF energies, surface energies for random alloys are much more challenging to compute, hence there are relatively few references available.One of the reasons is the presence of surface reconstructions, which can result in contracted or extended bond lengths near surfaces.Lattice spacing only gradually recovers to the bulk spacing, making it necessary to use relatively large cells for accurate calculations.In addition, finite size ordering effects on the surface must be avoided for random alloys.
According to Mak et al. [41], the {100} and {110} surfaces of MoNbTi exhibit surface energies of 2317 mJ/m 2 and 2171 mJ/m 2 , respectively.Moreover, the surrogate model of Hu et al. [27] has determined the surface energy for the {110} surface to be 2189 mJ/m 2 .Both results were obtained with the SQS approach.However, our methodology has produced lower values for these surface energies.Specifically, our calculations yield surface energies of 1989 mJ/m 2 and 1832 mJ/m 2 for the {100} and {110} surfaces, respectively.To understand the source of such a considerable discrepancy, we performed the con-vergence analysis with respect to the slab size (number of layers), which can be easily done with MTP.The analysis, whose results are presented in Fig. 9, clearly shows that the convergence of the surface energies can be rather slow.In the case of MoNbTi, up to 100 layers are required to produce the correct (converged) values.A more detailed examination of this behavior suggests that it can be traced back to a strong effect of statistical fluctuations of local atomic configurations on the relaxation contribution to the total energy (see SM S3 B for details).This emphasizes the importance of the analysis of supercell size convergence and also shows insufficiency of typical SQS sizes for accurate surface energies of random alloys.

B. Analysis of trade-offs
The VBA model is very efficient and accurately captures the strength and ductility parameters, enabling exhaustive scans across the entire concentration space.This allows us to track the evolution of material parameters as new alloy components are introduced and also quickly evaluate feasible regions with high resolution, which would otherwise imply significant computational costs.For example, by scanning over a dense mesh of concentrations for MoNbTi, MoNbTiTa, we can visualize the feasible regions of the two objectives (D and τ y ) as functions of the most influential material properties, with the results displayed in Fig. 10.For the planar defect energies, we use the dominating set of orientations ({110}/{112}), which determines the value of D for most of the Pareto-optimal compositions.
Examination of the feasible regions themselves already reveals some general trends.From the shapes of the regions one can infer that, generally, the CRSS correlates with the misfit parameter, δ, and the shear modulus, G, and only weakly with the USF energy, γ USF , while possessing practically no correlation with γ surf .Interestingly, despite the considerable correlation with the shear modulus, τ y exhibits no particular trends in terms of the Pugh's ratio, G/B.At the same time, the ductility index, D, has clear correlations with the surface and USF energies (in the form of the reciprocal relation) and a rather strong correlation with G/B, while the correlations with G and δ are relatively weak.An important conclusion to be drawn here is that, since τ y and D only weakly correlate with γ surf and δ, respectively, the latter two material parameters offer a possibility for an independent fine-tuning of the ductility and strength.
The trends inferred from the shapes of the feasible regions are, of course, not surprising.However, the analysis of the Pareto-optimal compositions provides us with a more refined picture.Firstly, one portion of the PFA for MoNbTiTa is almost completely overlapping with the PFA for MoNbTi.We will refer to this portion as MoNbTi(Ta).This portion of the PFA of MoNbTiTa corresponds to PFA1, which was previously discussed.The other visible portion of MoNbTiTa corresponds to PFA2.
The most prominent observation, when examining the Pareto fronts, is the consistent change in trends between the two portions of the PFA of MoNbTiTa.If one Pareto front shows a positive correlation, another one shows a negative correlation.For instance, although high surface energies are generally beneficial for ductility, D in MoNbTi(Ta) increases (i.e., the material becomes less ductile) at higher values of γ surf , mainly because γ USF (and G) increases faster.However, in the other portion of MoNbTiTa, this trend is reversed along the PFA, where the ductility improves (D decreases) thanks to a comparable rate of increase in both γ surf and γ USF .
Notably, the behavior of the Pugh's ratio, G/B, is much more non-trivial than what one could previously deduce from the plots of feasible regions.While it effectively distinguishes between brittle and ductile composition regions similar to D, the positive correlation with D is observed only along a more brittle part of the PFA for MoNbTi(Ta), while the more ductile part as well as the PFA of MoNbTiTa reveal a clear negative correlation with the ductility index.This suggests that blindly following the Pugh's criterion for improving ductility might result in very unfavorable compromises on the strength.
Finally, another way to use Fig. 10 is to rationalize the effects of the addition of Ta already visible in the Pareto-front representation in Fig. 3.The observed increased strength of optimal compositions with Ta (along the PFA) can be primarily attributed to Ta enhancing the misfit parameter, δ (see also Fig. 4a), and, to a certain extent, enhancing the shear modulus G.However, although high strength alloys are associated with larger values for both, G and δ, the maximum of these two quantities does not necessarily mean maximum strength, as is evident from the locations of the highest values of τ y at the Pareto fronts.Furthermore, the reduced ductility of Ta-containing alloys can be attributed to lower surface energies, which is not compensated by the modest decrease in USF energies.This also highlights how the Pugh's criterion can fail to predict the correct ductility trend when the alloy surface energy is ignored.
To check some of the above observations and to demonstrate the potential of the model for exploration of other systems, we extend our investigation of trends in D and τ y to a wider range of alloys comprising elements Nb, Ta, Mo, W, Ti, V, Zr, and Hf.Specifically, we consider the same alloy compositions as in Fig. 8 and apply the VBA model to compute all necessary properties for the main objectives, D and τ y , with the results shown in Fig. 11.We also add values from the Pareto set of MoNbTiTa (mind the difference in scales when comparing to Fig. 3), which puts this alloy system into a broader context.
A remarkable feature of Fig. 11 is an evident separation of points into two distinct clusters.These clusters can be categorized into specific alloy groups characterized by certain primary elements.One cluster, with points marked by red circles, consists of alloys rich in Ti, Zr, and Nb, which span a range of CRSS values from 200

C. Conclusion
In conclusion, we have demonstrated how a combination of CPA and actively-learned MTPs can solve the problem of predicting the CRSS and intrinsic ductility of refractory MPEAs over the whole composition space with ab initio accuracy, a problem that is currently not tractable with conventional ab initio-based methods.In particular, by integrating our ab initio automated workflow into a robust Bayesian multi-objective optimization framework we have applied the developed approach to identifying high-performance Pareto-optimal compositions.
Furthermore, our supercell size convergence tests enabled by MTP-powered molecular statics have revealed that sizes of typical supercell sizes used in conventional SQS calculations are often insufficient for accurate calculations of certain important materials properties, such as the surface energy, of multi-component alloys.Specifically, the surface energy of MoNbTi has been shown to exhibit a rather slow convergence with the supercell size, which leads to errors of 10-20 % for SQS supercells of about 100 atoms.
We have applied our multi-objective optimization framework to the four-component MoNbTiTa alloy system and exemplified its use by investigating how alloying with Ta affects the strength-vs-ductility tradeoff.Interestingly, we have found that, in this particular case, instead of advancing the Pareto front of the threecomponent system, the addition of Ta resulted in the emergence of a separate optimal region, extending the initial Pareto front of MoNbTi to higher strength values.Such a behavior represents a peculiar case of a disjoint Pareto sets in the composition space, which might have important implications for alloy design.
Furthermore, we have used the results from the workflow to validate a model based on the virtual bond approximation (VBA), which captures the bonding mechanism of transition metals.Relying on only a couple of element-specific parameters and universal functions, our VBA model has proved to provide estimates of alloy properties over the entire concentration space that are in very good agreement with direct CPA+MTP calculations.We have, then, used the VBA model to examine the trade-off between τ y and D in terms of relevant material parameters for a large class of refractory alloys.The results suggest that an improved trade-off in alloy design can generally be achieved if one fine-tunes the CRSS and ductility index by targeting the misfit parameter and the surface energy, respectively.This analysis also demonstrates how the VBA model can be a useful tool for examining trade-offs between competing properties in a large compositional space of refractory alloys.
Although our study focuses on strength and ductility, the proposed methodology can be readily applied to other mechanical properties, such as precipitation strengthening, which would involve interface energies or anti-phase boundary energies.Moreover, given the good accuracy of the proposed VBA model, it may also be used as a surrogate model in a multi-fidelity optimization framework, where ab initio evaluations of the objectives are triggered only rarely, allowing for efficient on-the-fly optimization of alloy properties.
Another prospective direction for development is the integration of recently developed multicomponent mag-netic MTPs [88] into the presented framework, opening the possibility to treat magnetic alloys with complex magnetic configurations.This can be also complemented by the learning approach, whereby the MLIPs are trained on the fly using fragments of the whole system [89], which could broaden the range of applicability of our framework to systems that are, in principle, barely tractable by conventional Kohn-Sham DFT, e.g., isolated dislocations.

A. Analytical model for intrinsic ductility
To predict the intrinsic ductility of a material, the Linear Elastic Fracture Mechanics (LEFM) based theory is employed.This theory, initially introduced by Rice and Thomson [57], has recently found application in the description of ductility of alloys and has been enhanced by incorporating atomistic insights [41,90].It aims to describe the behavior of ideal sharp cracks in elastic bodies and evaluates a material's propensity for failure due to crack propagation.This is done by balancing the stress intensity or strain energy release rate K Ic and K Ie of the governing processes.The competition between brittle cleavage and ductile dislocation emission at the crack tip determines the material's behavior, which is reflected in the stress intensity.Plastic deformation at the crack front results in a reduction of stress intensity, which ultimately can lead to an increase the total fracture toughness.We only consider in our investigation mode I loading (stress orthogonal to the local plane of the crack surface), which generally limits the ductility the most.Crack geometries that dominate ductility show cleavage on the {110} or {100} planes and dislocation emission on the {110} or {112} planes based on Ref. [41,42,91,92].We will refer to these fracture systems as {crack plane}/{emission plane}.The ratio of the K-factors for a given crack geometry gives the ductility index where intrinsic ductility is achieved for D < 1 at T = 0K.The stress intensity factor K Ic for cleavage fracture is given by and stress intensity factor K Ie for dislocation emission is given by where θ is the angle of the slip plane with respect to the crack front, ϕ is the angle of the Burgers vector b inclined to the slip direction, γ USF is the unstable stacking fault energy of the emission plane, γ surf is the surface energy of the crack plane and C is the linear elastic tensor.λ 22 is Stroh's anisotropic elasticity parameter [93].F 12 and o are geometrical and anisotropic parameters [90].

B. Analytical model for yield strength
The yield strength in bcc alloys with dominating solid solution strengthening can be evaluation with the model proposed by Maresca and Curtin [36,37].Within their model, the temperature and strain-rate-depended yield stress is characterized by two thermal activation models specific to certain temperature ranges.In the lowtemperature/high stress (τ y /τ y0 > 0.5) range, τ y (T ) is given by, whereas for higher temperatures/low stress (τ y /τ y0 > 0.5), In our calculations, we assume a ε of 5 • 10 −4 s −1 and ε0 of 10 4 s −1 .
In the above equations, τ y0 and ∆E b are, respectively, the zero-temperature yield stress and the activation barrier, given by the following expressions: ∆E b =2.00α where α = 1/12, B, µ and ν are averaged bulk, shear modulus and Poisson ratio, b is the average Burger's vector and the misfit parameter is δ The averaged linear elastic properties are obtained by: The misfit volumes ∆V i are obtained from a series of CPA calculations for systems with small deviations in the concentrations from the original alloy in the same way as in our earlier work [59].For each elemental component, four calculations are performed where the concentration of the respective element is changed by −δ, −δ/2, δ/2, and δ (δ = 0.005) from the original alloy.The concentration ratios among the remaining constituent elements are kept constant.The linear elastic constants were obtained from volume-conserving monoclinic and orthorhombic distortions following the computational details described in Ref. [94,95].To extract CRSS τ y from experimental data of poly-crystalline samples, we subtract an estimated Hall-Petch contribution from Ref. [96] and divide the result by the Taylor factor of 3.09.

C. First-principles evaluation of materials parameters
To determine the relevant materials parameters that affect ductility and strengthening, various ab initio methods based on density-functional theory are used.
Elastic constants, equilibrium volumes and concentration derivatives of equilibrium volumes are evaluated with the exact muffin-tin orbital (EMTO) [97] code (Lyngby version [98]) implementing a Green's function based DFT methodology combined with CPA to perform total energy calculations.Parameters for screened Coulomb interactions in the CPA are obtained by the use of the EMTO-based locally self-consistent Green's function technique [99,100].Accurate total energies are obtained within the full-charge density formalism [101].We evaluate finite-temperature equilibrium volumes following our previously suggested methodology [59] with LDA as the reference functional for the pressure correction calculation.The elemental reference volumes are obtained by extrapolation of cryogenic lattice constant measurements [102][103][104][105][106][107][108] to 0K.Zero point vibrations are included in the calculations by a quasi-harmonic Debye model U = 9/8N k b Θ D (V ) as described in Refs.[109,110].
The unstable stacking fault energies, γ USF , and surface energies, γ surf , are obtained using Moment Tensor Potentials (MTPs) [53,[111][112][113] which are actively-learned from DFT.The total energy of MTPs is given by a sum over all i atoms in the configuration where ξ α are linear parameters and B α (n i ) basis functions that depend on the neighbourhood n i within some cutoff radius around atom i.The scalar basis functions are obtained from contractions of the moment tensor descriptors where the per-atom radial functions f µ (|r ij |) are expanded in terms of Chebyshev polynomials.The total number of basis functions contained in E MTP is defined by the MTP level.For the hyperparameters of the MTP, we choose a level of 16 and a cutoff of 5.0 Å.The training sets are constructed using the algorithm proposed in [56].Below, we highlight only the most important points of the algorithm including our small modifications.For a more detailed description we refer to the original work of Hodapp and Shapeev [56].
In the first step, we randomly sample a set of training candidates containing 72/96-atom surface, 54-atom bulk, and 72/96-atom unstable stacking fault configurations on a regular concentration grid within a simplex.To explore more extrapolative configurations, we also introduce slight displacements of atoms and deformations of the cells.Choosing a grid with four points in each dimensions was found to be sufficient to predict energies over the entire composition space for the initial training set (for details, see [56]).In total, our initial training set contains 146 and 377 configurations.Next, we let active learning choose the most distinct configurations from the set of training candidates (each candidate set contains around 10000 configurations), perform static DFT calculations on them and augment the interatomic potential by retraining it with these data.The search for distinct configurations is conducted using a progressively decreasing extrapolation grade, starting from 1000 and gradually reducing to 100, 10, and finally 2. We iterate multiple times, with always newly generated training candidate sets and stop when no new training configurations emerge anymore that exceeds the extrapolation grade.This procedure resulted in a training set containing 219 and 631 configurations for MoNbTi and MoNbTiTa, respectively.Finally, we run molecular dynamics simulations at 50K, 150K, and 300K, on each of the configurations from our initial training set.During each of the simulations, active learning is used to find new potentially extrapolative configurations.We remark here that running MD was crucial for the robustness of our potentials.Training only on relaxations as done in [56] has been found to be insufficient to ensure stability of large-scale configurations during the lattice relaxation.Our final training set contains 381 and 720 configurations for MoNbTi and MoNbTiTa, respectively.The errors on the final training set for total energies are found to be 6 and 19 meV/atom, and for forces -0.099 and 0.104 eV/ Åfor the MoNbTi and MoNbTiTa potentials, respectively.A detailed performance check can be found in SM S3 C.
The final properties (unstable stacking fault and surface energies) are then calculated on a cell with about 300000-500000 atoms to include far-field effects and to ensure proper statistics.We provide a detailed convergence study for the cell size in SM S3 D.
The DFT computations for the MTP training are performed with the Vienna Ab-initio Package (VASP) [114][115][116] using the PBE functional [117] with an energy cutoff of 400 eV.We use a Monkhorst-Pack k-point mesh with a spacing of 0.15 Å−1 .Moreover, we use Gaussian smearing with a smearing width of 0.08 eV to avoid potential problems of negative occupations.The standard pseudopotentials (PP) were used for Mo, Ti, and Ta [118,119].For Nb, PP including semicore p-states as valence states was used.

D. Multi-objective Bayesian-optimization
In order to effectively find the Pareto front of our multi-objective Bayesian optimization problem, we have adapted a methodology developed by Galuzio et al. [120].This method involves generating a Pareto front approximation (PFA) from an easy-to-evaluate surrogate model at each search step, and selecting the next evaluation point from this PFA before updating the surrogate model.This iterative process enables the PFA to gradually approach the true Pareto front.We can also narrow the search space by excluding compositions that are thermodynamically unstable or possess unfavorable properties.
To navigate through the search space, we employ the following sequential design strategies: 1. We begin the search process by building a smooth surrogate model.To that end, we are using Gaussian processes regression (GPR) based on the observations available from previous rounds.For the initial round, we sample a coarse regular grid in the search space.The choice of the covariance function K(x, x ′ ) is determined using cross-validation on grid test sets.We have selected a sum of white noise kernels and radial basis function kernels as the best: where σ, C, λ are hyper-parameters estimated by maximizing the marginal likelihood on the training data.
The surrogate models for our objectives, namely the ductility index D and the yield strength τ y , are then used to make predictions at unobserved areas of the search space and quantify the uncertainty around them.This information is crucial for guiding the search towards promising areas in the search space.
2. Next, we use the non-dominated sorting genetic algorithm II (NSGAII) [58] implemented within the DEAP library [121] in order to identify 100 PFA points Ω PFA using our surrogate model.A penalty function is used to handle constraints, such as limiting the search space and keeping it within a certain concentration range.
3. Select the next search points x n+1 from the current PFA, x n+1 = x * tnext ∈ Ω PFA , that are farthermost from all previous sampling points according to where {•} corresponds to the smallest distance from the i th point in Ω P F A to all other points found by the algorithm.The subscript "f" indicates the distance cal-culated in objective space, while the subscript "x" indicates the distance in the search space.The distances are standardized by µ {•} , representing the average, and σ {•} , representing the standard deviation, of the respective distance measure.We choose the weighting parameter q as 0.8 to slightly favour objective space exploration.
4. Check the stopping criterion: If the maximum number of iterations has been reached or only minimal Pareto-improvement is achieved, stop the algorithm and output the Pareto front and the Pareto set.The stationarity condition on the Pareto front measures Paretoimprovement by determining how much the front changes with more points.If there is little change after several iterations, the loop terminates.Otherwise, go back to step 3 and continue the search.The final Pareto set can be discontinuous and portion of the Pareto set can be located in different regions of the concentration space.To identify separate portions, we use multivariate Gaussian mixture models to perform automatic spatial clustering and label the corresponding parts at the Pareto front.For MoNbTi, 36 points are required, whereas for MoNbTiTa, 143 points are needed, with 72 for set 1 and 71 for set 2.

E. Dimensionality reduction and visualization
We employed a t-Distributed Stochastic Neighbor Embedding (t-SNE), a dimensionality reduction technique, to visualize the 4-dimensional concentration space in two dimensions.This projection allows us to visualize our objectives and gain insights into the overall landscape.The biggest advantage of t-SNE is that it can capture much of the high-dimensional data's local structure while also revealing global structure.Furthermore, it is a non-linear method that can also capture the structure of complex manifolds.However, t-SNE provides only an approximate preservation of local structure and distances, which leads to the emergence of an observed patchy pattern in our case.The corner points of the projection correspond to the elemental compounds, the edges represent binary compositions, and the center represents the equimolar compound.The distance from the corner points in the projection is approximately proportional to the effective composition.The t-SNE projections displayed in Fig. 3 in the main text are performed from the entire concentration space, with points evaluated using the Gaussian process surrogate models for the ductility index and the strength.

F. Virtual bond approximation model
Within the VBA, to second order in inter-component interactions, the energy of an alloy system can be ex-pressed as follows: where c i are atomic fractions of components, e (k) are the virtual bond functions of the average d-valences, N (k) ≡ {N i , N ij }, and of bandwidth parameters w (k) ≡ {w i , w ij }.Here N i and w i refer to a respective elemental value, while second-order parameters are defined as N ij = (N i + N j )/2 and w ij = √ w i w j , where the latter expression is inspired by the tight-binding description of transition-metal alloys [122].λ represents all other parameters, such as temperature, volume, lattice vectors, etc.The bandwidth parameters were taken from Ref. [123] and adjusted by a small correction depending on the row of an element.Alternatively, one can fit these parameters to elemental compounds and a few binary alloys but usually the relative values obtained this way are similar to those from the above reference.Note also that the expression E({c i }, λ) can represent either the Gibbs free energy, G({c i }, P, T ), or the Helmholtz free energy, F ({c i }, V, T ), where additional parameters (e.g., lattice vectors) are implied.Considering derivatives of the thermodynamic potentials with respect to various parameters, we can get VBA representations of any property that can be derived from the free energy.
From the Gibbs free energy, we can the derive an expression for the volume, V ({c i }) = ∂G({c i }, P )/∂P | P =0 , omitting the temperature dependence for simplicity.Similarly, we can obtain an expression for the bulk modulus and general elastic (isothermal) constants by , where ϵ α is the strain in Voigt notation.In the context of surface and USF energies, the molar energy difference, ∆E = E defect − E bulk , between the defective and bulk structures is used complemented by the area normal to the planar fault deduced from the predicted volume.The final expression always assumes the form: where ξ({c i }) is one of B, C ′ , C 44 , V , or ∆E, and ξ (k) represent the corresponding virtual bond functions.To obtain the unknown ξ (k) , we perform a series of calculations of the associated quantity for several systems (elemental compounds and selected binary, three-and fourcomponent alloys).The unknown values of virtual bond functions at integer (for ξ (1) ) or half-integer values (for ξ (2) ) are then obtained by means of linear regression from the above equation, since the values of c i and w i are known for each composition.The discrete function values obtained in this manner for our refractory alloys behave very well and can be represented as values of smooth functions of valences.In fact, we use second and thirdorder polynomials to accurately represent, respectively, the first order virtual bond function ξ (1) and the second order virtual bond function ξ (2) as a way of conveniently obtaining a value for a given N .This implies that virtual bond functions can be used to consistently characterize the system.An example of the virtual bond function can be found in Supplementary S4.
The material parameters for determining the corresponding virtual bond functions have been obtained using CPA calculations.For planar defects we first calculate unrelaxed surface and USF energies with CPA for elemental compounds and binary alloys and parameterize the VBA model, as described above.We have checked that the VBA model predicts the CPA values very well (see Supplementary S4).To account for relaxation effects, we simply adjust the virtual bond functions to align with the relaxed surface and USF energies of individual elemental compounds.More precisely, we scale the virtual bond functions by a factor that follows a linear relationship with the d-valences.Such an approach works because fitting to CPA values already includes most of the chemical interactions and the strongest relaxation effects have geometrical nature, which can be taken into account by the rescaling procedure of the universal virbond functions.
Supplemental Materials: Ab initio framework for deciphering trade-off relationships in multi-component alloys

S1. VALIDATION OF PHASE STABILITY OF SOLUTE SOLUTION
To assess the stability and formation of solute solutions and the occurrence of competing intermetallic phases, various approaches can be employed.
One approach involves CALPHAD utilizing thermodynamic databases, albeit constrained by the availability and also accessibility of such databases.Alternatively, there are also fully ab initio-based approaches to evaluate phase stability, although they tend to be computationally too expensive and complex to be utilized inside an automated workflow.Accurate and reliable prediction of phase transformations can be a non-trivial task even for binary systems.
A more straightforward strategy involves assessing simpler quantities such as the formation energies of solid solutions and ordered phases, complemented by exper-imental observations to establish criteria for estimating phase stability.
We employ the formation energies criterion from Ref. [124] in our work to evaluate operational stability of our alloy and to set a priori boundaries of the design space.This criterion relies on the observation that multiple phases form when E f > 70 meV and metastable metallic glasses form or ordered phases appear when E f < −150 meV [125].In our case, the formation energies of the entire system fall within these boundaries, as illustrated in Fig. S1 and Fig. S2.
Furthermore, we also take into account the ordering energies for binaries on or near the convex hull of formation energies.The ordered phases were obtained from supplementary materials of Zheng et al. [61] and materialsproject.org.The ordering energies are determined by subtracting the formation energy of the ordered intermetallic phase from the formation energy of the solid solution, In Fig. S1, it is evident that only the MoNb binary side exhibits a significant ordering energy, peaking at around 80 meV.In the four-component system, slightly higher ordering energies are observed at the TaMo binary side, reaching 120 meV.However, this corresponds only to a small region close to the boundary of the concentration space, which is irrelevant for the optimization.Moreover, in a multicomponent alloy, nucleation of an ordered phase consisting of two elements can be hindered by interactions between these two elements and other components.
Most importantly, experimental evidence strongly supports the assumption that MoNbTi and MoNbTiTa form solid solutions.Mo, Nb, and Ta are fully miscible at typical annealing temperatures and readily form a solid solution.Various experimental studies have demonstrated that the addition of Ti to the MoNbTa system also leads to the formation of a BCC solid solution [126][127][128].−0.10 −0.05 0.00 0.05 0.10

S2. VALIDATION OF MECHANICAL PARAMETERS A. Validation of Ductility calculations
Experimental verification of the ductility model for a variety of refractory MPEAs has been performed in Ref. [41].Following the same approach, we compute the ductility index, D, as a minimum of two values, D 100/110 and D 110/112 , corresponding to two surface orientations, {100} and {110}, and two stacking fault orientations, {110} and {112}, respectively.USF, γ U SF , and the surface energy, γ surf , is obtained from molecular statics of fully relaxed structures using MTPs.Unlike the SQS approach, which requires full structure relaxation calculations for each property, crystallographic orientation, and composition in a cell with typically 800-1000 atoms for a three-component alloy, our active-learning approach generates 381 and 720 static configurations with 54-96 atoms/configuration to fit an accurate potential in the entire composition space of MoNbTi and MoNbTiTa, respectively.
Using the USF and surface energies, along with the linear elastic parameters from CPA, we evaluate parameters D 100/110 and D 110/112 and ultimately the ductility index, D. Our values of the ductility parameters, D 100/110 = 1.53 and D 110/112 = 1.35, in the equimolar MoNbTi, corresponding to x = 0 in Fig. S3, are somewhat different then the ones (1.29 and 1.59, respectively) reported in Mak et al. [41].This can be attributed to supercell size convergence issues associated with SQS surface energy calculations.In Fig. S3, we examine D for the MoNbTi alloy as a function of varying concentrations of each component A ∈ {Mo, Nb, Ti} in the two fracture systems, 100/110 and 110/112.The concentration of A is varied, keeping the proportion between the rest of the alloy equal, i.e., the general formula of a modified alloy is Notably, the 110/112 fracture system is the dominating one across all cases, as it has the lowest D index.Increasing Nb content enhances ductility, while Mo and Ti additions show no positive effect.Nb's softening properties are experimentally documented in Refs.[129][130][131].

B. Validation of solid solution strengthening calculations
Maresca and Curtin [37] suggested that in concentrated alloys, especially at high temperatures, solid solution strengthening can be described based on the alloy's average volume, linear elastic moduli, and misfit volumes of components.In particular, the determination of linear elastic constants and concentration derivatives of equilibrium volumes becomes challenging when employing supercell methods such as Special Quasirandom Structures (SQS).The presence of broken symmetry in SQS intro-duces direction-induced errors in elastic constants that are challenging to alleviate [47].Additionally, obtaining concentration derivatives poses difficulties, as SQS is not effective for arbitrary concentrations.All of these parameters can be obtained robustly within CPA [59,132,133] because local atomic relaxations have a minor effect.
Given the computational efficiency of CPA it eliminates the need for extrapolative models for misfit volumes, distinguishing it from prior studies [16,32].While misfit volumes often exhibit a linear dependence on concentration in refractory alloys, our research identifies specific alloys where this relationship deviates notably.Fig. S4 illustrates the misfit volumes and bulk modulus that were directly calculated alongside those determined using ROM for variations of Ta within the MoTaW alloy.Similarly, misfit volumes and bulk modulus shows quite significant deviations form the linear trends.Especially around 0.28 the bulk modulus abruptly changes, which can be attributed to topological transitions.In Fig. 4b, we find good agreement between our CPA results and those obtained from an embedded-atom potential for Aatom in larger cells, as detailed in Refs.[72][73][74].
To verify the final accuracy of the strengthening calculations, we compared the accuracy of our CRSS values calculated at 300 K to room-temperature (RT) experimental data from tensile tests [70,71], as well as compression tests reported in Ref. [65].The poly-crystalline samples used in the experiments were reported to be singlephase.To extract the CRSS from these tests, we subtract the estimated Hall-Petch contribution of 20 MPa [96] and divide the result by the Taylor factor of 3.09.In comparison to the experimentally measured total yield strengths, the strengthening contribution of the grain boundary is small.
However, for NbTaTi, the deviations are quite large.Our theoretical τ y value for NbTaTi is approximately 88 MPa, which is about half of the experimental value obtained from tensile tests, which yields a CRSS of around 168 MPa.To elucidate the nature of this discrepancy, we have performed a comprehensive analysis of solid solution strengthening in five additional 3-and 4-component MPEAs: NbTiZr, MoNbTi, CrMoTaTi, MoNbTaW and MoTaW.Among the alloys tested, NbTiZr, MoNbTi, and MoNbTaW exhibit similar τ y values of around 300 MPa in experiments, whereas our theoretical approach yields 293, 195, and 234 MPa, respectively.The 4-component CrMoTaTi alloy, which exhibits a yield strength of 621 MPa in experiments, displays close with our theoretical prediction of 577 MPa.
By comparing Fig. 4a with Fig. 4c, one can see that the theoretical predictions aligns closely with experimental data when the average misfit δ is high, as exemplified by NbTiZr and CrMoTaTi.All other alloys have smaller misfit deltas and their theoretical predictions consistently underestimate strength.Furthermore, it seems that the relative error is reduce for stiffer alloys with larger elastic moduli.
Calculating CRSS in bcc crystals is challenging due to multiple active glide mechanisms.According to Maresca and Curtin [36,37], the solid solution strength of our multi-component refractory metal alloys is to large extend influenced by interactions between edge dislocations and solutes, with a lesser contribution from solute-screw dislocation interactions.Building on that strengthening theory, Baruffi et al. [134] discovered that as the misfit parameter δ or elastic moduli increase, the dominance of edge-controlled contributions becomes more pronounced, eventually determining the strength completely.Within the used analytic strengthening model, we, however, only consider the edge-controlled strengthening.The contributions arising from screw dislocation are omitted, because their contributions are very challenging to calculate.The negative example of NbTaTi has a relatively small misfit parameter δ and is also soft.As a result, we underestimate the strength because the screw-controlled contribution is significant and can actually not be neglected.In contrast, NbTiZr and CrMoTaTi are wellpredicted because the screw contribution can be effectively disregarded.We want to emphasize, however, that as per Maresca's findings [36,37], the edge-controlled contribution tends to become increasingly dominant at elevated temperatures.Given that refractory alloys are primarily intended for high-temperature applications, this observation strengthens the validity of using a model for strengthening that is solely based on edge contributions for all alloys.Even though we only consider a portion of the strengthening contribution, the findings presented in Fig. 4c demonstrate a high level of consistency between the experimentally determined CRSS and our calculated values.Overall, we tend to underestimate the strength of all the alloys tested, as expected, because we ignore the screw contributions.This also implies that, we can use the values for a conservative estimate of strength.In general, the relative changes in strength that occur during alloying are accurately reproduced.We can distinguish a weaker and a stronger alloy.For example, hardness screening of Nb in MoNbTaW revealed only a modest strength increase of about 10% [135] from equimolar MoNbTaW to MoTaW (no reliable yield strength measurements are available for MoTaW).Our methodology predicts a similar magnitude of increase, from 234 MPa to 257 MPa.

S3. VALIDATION OF MATERIAL PARAMETERS A. Convergence of local relaxation on surface and stacking fault energy
Interatomic potentials have the advantage of allowing for larger cell sizes compared to SQS methods with.Using large cells is especially important when studying farreaching elastic effects that cannot be accurately captured in small cells or need some complex compensation procedures.In the case of surface energies, the defect structure was created by introducing a vacuum of 20 Å.We checked that fixing atomic positions of certain deep internal layers had no effect.For both surface orientation, we create orthogonal cells.In case of USF, we created a defect structure with two stacking faults in cell by simply shifting half of the atoms in 1  4 a [111].We keep the cell orthogonal.We only allow for atomic relaxation normal to the fault plant and do not change the out-of-plane lattice dimension.By using large enough structures with more than 100 layers, we can account for inelastic displacement associated with the fault.This is in contrast to [41,90].
In order to assess the necessary cell size for the calculation of surface energies, we conducted convergence tests.Fig. 9 displays the surface energies for the 100 orientation as a function of the number of layers perpendicular to the fault.Convergence varied by alloy, and MoNbTi exhibited the slowest convergence.Approximately 100 layers were sufficient for convergence for the more sparselypacked 100 layer in the case of MoNbTi.Only about 50 110 layers are required to achieve convergence.The other alloys, MoNb and MoNbTa, require 8-10 layers to achieve convergence.
On the other hand, the necessary cell size for unstable stacking fault converges much faster, even after only a few layers.We want to stress that we don't perform optimization of the neighboring cluster vector.Our big cell will contain statistical clustering of atoms of a specific kind.We assume that by randomly distributing the atoms in the cell in a sufficiently large cell size, we can generate a truely random alloy.For comparison, we also performed reference calculations with 4x4x10 SQS cells with 160 atoms with cluster vectors closely resembling those of truly random alloys.Except for MoNbTi, the energies of SQS cells are numerically equivalent to those of large cells.Numerical errors related to the interatomic potential can be assumed to be below 5% [56].one alloy component in relation to the equimolar composition, while keeping the ratios of the remaining components constant.Neglecting the effects of local relaxation on USF γ USF yields a consistent increase of about 25 % for both orientations.The overall trends of γ USF remain unchanged upon changes in composition.For surface energies γ surf , the impact of relaxation is significant and varies depending on the element in question.Without considering relaxation effects, surface energies are largely unaffected by the addition of Nb.However, when relaxation effects are taken into account, the addition of Nb causes a significant increase in γ surf .Similarly, changes in Ti concentration lead to more pronounced effects on γ surf .On the other hand, for Mo, there is an almost constant shift between the relaxed and unrelaxed cases.
To show the origin of this behaviour, we used the Polyhedral Template Matching (PTM) [136] method to identify the local crystalline structure in order to relate this relaxation contribution to local structural changes.It also provides for each atom in the structure a RMSD value, which is a measure of the spatial deviation from the ideal local structure.The distribution of the RMSD value for each atomic species in the structure is Gaussian-like.Especially Ti showed higher local distortions compared to all other components in all alloys investigated.We also found the mean and standard deviation of the RMSD correlate to the energy difference between unrelaxed and relaxed surface energies.The higher the standard deviation of these RMSD values, as seen in Fig. S7, the more significant the relaxation contribution will be.The overall effect of this local structural changes, seems to be less relevant for USF.

C. Performance metric for MTPs
The surface and unstable stacking fault energies are relatively small and sensitive quantities.Even minor errors can have significantly impact on the final results.Therefore, simply providing errors on forces and total energies calculated with MTP for target structures may be inadequate to demonstrate the accuracy of the potential for a particular property.To verify accuracy of γ surf and γ USF , we generate random structures of surface and unstable stacking fault configurations, along with their respective bulk reference configurations.Atomic species are placed randomly on lattice sites for various alloy compositions and the atomic positions and the cell volumes are slightly randomly perturbed.Atomic relaxation is then performed using our MTPs.The reference energies are then computed on the relaxed structures using DFT.Across the test set, mean absolute difference between the MTP and DFT forces is found to be 0.0402 eV/ Å.Additionally, we provide a parity plot in Fig. S8 to compare the MTP and DFT surface and unstable stacking fault energies.Overall, there is a high degree of correlation between MTP and DFT proving the accuracy of our MTPs.

D. Influence of the elemental distribution
In order to estimate the influence of the elemental distribution, we investigate the distributions of surface γ surf and unstable stacking fault energies γ USF using atomic configurations with different distributions of atomic species.To that end, we generated several hundred configurations for MoNbTi with fully random atomic distribution on the lattice sites, and using SQS.For SQS, we optimized the pair correlations of the first four shells, and the first triplet cluster correlations.We perform the optimization of the atomic configurations on the bulk structures and neglect 2-dimensional correlation at the defects.Fig. S9 illustrates the expected narrowing of energy spread with increasing cell sizes.The standard deviation in the energies is smaller for the unstable stacking faults than for surface energies.Furthermore, the energy distribution of the SQS configuration for the stacking faults show a similar standard deviation to that of substantial bigger cell of around 10000 atoms with a purely random atom distribution.This result is in agreement with the distributions reported in the supplementary material of [61], who also used MTPs.

E. Error propagation
Two primary error sources arise: first, from representing a truly random alloy within a finite cell, accounting for variations in atomic configurations; and second, from approximations made by the potential.
From Fig. S9 we deduce that the distributions for the surface and unstable stacking fault energies due to different atomic configurations are close to Gaussian distribution.Hence, we can perform an analytic Gaussian error propagation to estimate the final error on the D parameter, Let us denote the errors as δD, δγ USF , and δγ surf .Since χ is computed from CPA values, it is not subjected to statistical errors and can be treated as a constant.Assuming statistical independence of δγ USF and δγ surf , we have with the partial derivatives, So, the error in D can be written as Using the expectation values and standard deviations for the largest cell size, we obtain a value for D of 1.35 ± 0.004.Hence, we conclude that taking a large enough cell size is sufficient to practically eliminate any influence of the elemental distribution on D.
Obtaining the second type of error, originating from the potential, is more challenging due to limitations imposed by the size of the test set.If one were to consider the Root Mean Square Error (RMSE) values from Fig. S8 as error estimates and assume complete independence of the errors, a straightforward calculation would yield an error of δD = 0.02.This value is approximately five times larger than the influence of the atomic contribution.However, it's essential to acknowledge that the errors stemming from the potentials of γ USF and γ surf are, in fact, correlated.
To account for this correlation we add the covariance between γ USF and γ surf in error propagation: (δD) Incorporating the covariance, we find δD = 0.015.It is important to note that this error represents a conservative upper bound.For instance, it is clear from Fig. S8 that the largest contribution to the average error in γ surf comes from the region of very high surface energies.Since alloys in this region are also characterized by large γ USF they are expected to be very brittle and hence far from the Pareto front.Consequently, the overall error is small enough to accurately calculate the D-parameter and its concentration dependency.Moreover, smaller deviations encountered during the Pareto front search are further smoothed out through the appropriate choice of the Kernel for the Gaussian process regression and appropriate sampling multiple configuration.

S4. VALIDATION OF VBA MODELS
A. Validation of the virtual bond functions Fig. S10 shows virtual bond functions, v (1) and v (2) , fitted for the volumes.They are smooth function of the valences and are fitted with second order polynomials.The misfit volume is directly derived from the volume expansion as follows: The volume expansion is primarily influenced by the first-order parameter function, v (1) .The second-order parameter function v (2) can lead to non-trivial contributions to the equilibrium volume and its concentration derivatives in multi-component alloys.These functions' absolute values are determined by the scaling of the rowdependent prefactors and the band-width factor itself.The row-dependent prefactor for the bandwidth parameters is optimized for the training data of each individual property.The virtual bond function are not necessarily linear for other properties.v (1)  v (2)   FIG.S10.Virtual bond functions v (1) and v (2) fitted for the volumes.that shows significant deviations, especially for decreasing fraction of Mo.This is especially true for alloy systems such as MoNbTi, where relative errors between the VBA model and the rule-of-mixture model are noticeable already for the equimolar system.However, in the case of the elastic constant, specifically C 44 , noticeable deviations become apparent when employing the rule-of-mixture approach (as illustrated in Fig. 5a).While the rule-of-mixture does manage to capture the general trend, the discrepancy increase significantly, particularly for the 4-component systems and MoTaW, where it can reach 40 GPa.In contrast, the VBA model showcases good agreement with direct calculations for the majority of the alloys, exhibiting only minor discrepancies for NbTiZr and NbTaTi.
So far we have showed results for quantities that can be directly obtained from CPA calculations.However, we cannot apply the same approach to the surface and USF energies because of the appreciable effect of atomic relaxations.To avoid the need of heavy SQS calculations or of training MTPs for several test alloys we take a simpler route.First, we calculate unrelaxed surface and USF energies with CPA for elemental compounds and binary alloys and parameterize the VBA model, as described above.We have checked that the VBA model predicts the CPA values very well.To incorporate relaxation effects, we simply adjust the virtual bond functions to align with relaxed surface and USF energies of individual elemental compounds.More precisely, we scale the virtual bond functions by a factor that follows a linear relationship with the d-valences.Such an approach works because fitting to CPA values already includes most of the chemical interactions and the strongest relaxation effects have geometrical nature, which can be taken into account by the rescaling procedure of the universal virtual bond functions.
In Fig. S12, we compare directly computed surface energies from CPA (without relaxation) and PAW (with relaxations) to our model's estimates using ROM and virtual bond approximations.The energies of surfaces and the unstable stacking fault are effectively characterized by VBA model, which outperforms the simpler ROM.Notably, the VBA model significantly enhances the precision of surface energy predictions, underscoring its efficacy in capturing intricate energy behaviors.The corresponding results for the {110} surface energy varying with the concentration of individual components are shown in Fig. S13, from which one can see that the simple rescaling is sufficient to take into account the concentration dependence of relaxation effects.
In Fig. S14, we extend our comparative analysis to the elastic constants, following a similar approach.Notably, when examining the elastic constant C 44 , it becomes evident that the VBA demonstrates a substantial enhancement compared to the traditional ROM.This improvement underscores the efficacy and accuracy of VBA in capturing essential material properties beyond what can be achieved through conventional methods.

FIG. 1 .
FIG. 1. Flowchart of the materials design loop using Bayesian multi-objective optimization: (a) Advanced ab initio techniques such as the coherent potential approximation (CPA) and actively-learned moment tensor potentials (MTPs) to calculate the strength, τy, and ductility, D, objectives.Quantities obtained from ab initio are: Shear, C44, C ′ , and bulk, B, moduli, specific volume (per atom), V , and its concentration derivatives, ∂V /∂c (n) , all evaluated at room temperature; the surface energy, γ surf , and the USF energy, γUSF.(b) Employing Gaussian process regression, we construct a surrogate model, and a genetic algorithm is applied to the search of Pareto-optimal compositions.(c) The sampling points are chosen from the Pareto front approximation until convergence is reached.(d) Pareto-optimal alloys are the key outcome of the process.

FIG. 2 .
FIG. 2. Optimization results for the MoNbTi system: (a) Pareto set approximation (PSA) alongside with sampling points in the search space.(b) Pareto front approximation (PFA), sampling points and feasible region in objective space.PFA was gradually color-coded from yellow, which corresponds to high strength, to pink, which corresponds to high ductility.PSA uses the same color coding.

FIG. 3 .
FIG. 3. Optimization results for the MoNbTiTa system: (a) Visualization of the PFA and sampling points.The Pareto set comprises two distinct regions in the concentration space derived from automated clustering.These regions define the corresponding sections (PFA1) and (PFA2) on the Pareto front.The sampling points are color-and marker-coded based on their proximity to the respective PFA regions.Open symbols mark specific alloys mechanically tested in Refs.[65, 66].t-SNE projections of (b) ductility index D and (c) strengthening τy for the whole search space.Grey areas are beyond the search space.Sets 1 and 2 mark sampling points that are closer in the search space to the points of PFA1 and PFA2, respectively

FIG. 6 .
FIG. 6. Correlation between directly calculated D and τy and their counterparts DVBA and τy VBA derived from the virtual bond approximation.

FIG. 7 .
FIG. 7. Parity plot of predicted and experimentally obtained CRSS, τy, for various alloys, including alloys containing Hf and V.

67 FIG. 8 .
FIG.8.Correlation plot of predicted D and fracture strain for various alloys from Ref[23].Alloys containing Hf or V are marked with orange and brown dots, respectively.Alloys containing both Hf and V are marked with pink dots.

FIG. 9 .
FIG. 9. Convergence of surface energies with respect to the number of layers for pure Mo, equimolar MoNb, MoNbTi, and MoNbTa, along with corresponding results for an SQS cell sized 4x4x10.All values are obtained with MTPs.

FIG. 10 .
FIG. 10.Feasible regions (light blue for MoNbTi, light red for MoNbTiTa) and correlations between materials parameters (surface energy, γ {110} surf , USF energy, γ {112} USF , (c) shear modulus, G, misfit parameter, δ, and Pugh's ratio, G/B, where B is the bulk modulus) and objectives (the ductility index, D, and CRSS, τy) evaluated along the Pareto fronts (marked with blue-green for MoNbTi and with red-yellow for MoNbTiTa).MoNbTiTa's PFA has two segments: PFA2 is represented by the predominantly yellow tail, PFA1 overlaps with the PFA of MoNbTi.

FIG. 11 .
FIG. 11.Ductility index, D, and CRSS, τy evaluated for various alloys using the VBA model.Red circles indicate Ti, Zr, and Nb-rich alloys, while blue circles denote Ta and W-rich alloys, with the fill color (from dark blue to yellow) corresponding to the measured fracture strain ([23]).The Pareto front of MoNbTiTa is plotted in magenta points.
FIG. S1.Ternary plot showing a heatmap of the formation energies of the MoNbTi solid solution.Dots denote binary compositions corresponding to intermetallic compounds on the convex hull from Ref. Zheng et al. [61], with the fill color indicating the respective ordering energies.

Mo 100/ 110 Nb 100 / 110 Ti 100 / 110 Mo 110 / 112 Nb 110 / 112 Ti 110 / 112 FIG
FIG. S3.Comparison of the ductility index D for the different fracture systems 100/110 and 110/112 (denoted by subscript) upon changing the concentration of one component in the equimolar MoNbTi alloy according to Ax[MoNbTi]1−x, where x is the molar fraction added to the equimolar alloy.
FIG. S4.Comparison of directly calculated and valuesthe rule-of-mixture for (a) misfit volumes and (b) bulk modulus of a MoTaW alloy with respect to Ta variation.

B. 2 ]
FIG. S5.Comparison stacking fault energy {110} and {112} with and without local relaxations upon changing the concentration of one component in the equimolar MoNbTi alloy according to x Y → MoNbTi, where x is the molar fraction of Y added to the equimolar alloy.

2 ]
FIG. S6.Comparison surface energy {110} and {100} with and without local relaxations upon changing the concentration of one component in the equimolar MoNbTi alloy according to x Y → MoNbTi, where x is the molar fraction of Y added to the equimolar alloy.
FIG. S7.Energy difference of unrelaxed and relaxed surface energies and standard deviation of RMSD value from polyhedra template matching for different compositions of MoNbTi.

2 ]
FIG. S8.Parity plot of surface and unstable fault energies obtained from MTP and directly from VASP DFT code.

B
Fig. S11 displays a comparison between direct and VBA calculations of misfit volumes, where the concentration of Mo in MoNbTi is varied away from the equimolar value.The figure illustrates that the VBA model is capable of capturing the trends in the concentration dependence very well, in contrast to the simple rule-of-mixture FIG. S11.Comparison of misfit volumes in MoNbTi withvarying concentrations of Mo obtaind from direct calculation (⋆, ×, +), rule-of-mixture (ROM) and virtual bond approximation (VBA).
FIG. S12.Comparison of surface energies obtained from Vegard's law (V), from virtual bond approximation (M) and calculated directly with CPA and PAW (C).PAW is relaxed.CPA unrelaxed.Virtual bond function are not rescaled.
Histograms of surface γ surf and unstable stacking fault energies γUSF for purely random configurations of different number of atoms and SQS generated cells obtained from MTPs of a near equimolar MoNbTi alloy.

TABLE S1 .
[23]arison of the fracture strain and D parameter calculated from the VBA model.Fracture strain values are from Ref.[23].