Impact of local lattice relaxations on phase stability and chemical ordering in bcc NbMoTaW high-entropy alloys explored by ab initio based machine-learning potentials

Recently, high-entropy alloys (HEAs) attracted wide attention due to their extraordinary materials properties. A main challenge in identifying new HEAs is the lack of efficient approaches for exploring their huge compositional space. Ab initio calculations have emerged as a powerful approach that complements experiment. However, for multicomponent alloys existing approaches suffer from the chemical complexity involved. In this work, we implement a new method for studying HEAs computationally. Our approach is based on the application of machine-learning potentials based on ab initio data in combination with Monte Carlo simulations. The high efficiency and performance of the approach are demonstrated on the prototype bcc NbMoTaW HEA. The approach is employed to study phase stability, phase transitions, and chemical short-range order. The importance of including local relaxation effects is revealed: they significantly stabilize single-phase formation of bcc NbMoTaW down to room temperature. Finally, a so-far unknown mechanism that drives chemical order due to atomic relaxation at ambient temperatures is discovered.

approaches consistently showed such an incipient B2 ordering when the alloy is cooled down from the solid solution. 11 A limitation of such perturbational approaches is the limited inclusion of local relaxation effects. Indeed, explicit supercell calculations of a B2 ordered state as well as a disordered solid solution revealed that the order-disorder transition temperature can be significantly lower if local lattice relaxations are included in the computations. 12 Recently, Wang et al. extended this approach and considered in total 178 ordered supercell configurations to study the phase stability of bcc NbMoTaW. 13 The B2 ordered structure has been included in the pool of structures and it was found to be stable below 600 K. The limitations of such approaches are, however, that the considered structures must be anticipated, i.e., they must be included in the data pool a priori. Moreover, ideal mixing of elements is assumed not accounting for possible short-range order effects.
The CE can include implicitly local lattice relaxations and can be combined with Monte Carlo simulations to account for chemical short-range order. CE, however, becomes computationally demanding when studying systems with a large number of chemical elements 14 due to increased combinatorial complexity of possible interatomic interactions. Therefore, in practice, for multicomponent alloys with more than three elements the amount of interactions is typically limited to a few interactions. For example, based on nearest-neighbor pair interactions, Huhn and Widom 15 also found a B2 ordering at ambient temperatures. Due to these limitations it is, however, not clear whether such an ordering at ambient temperatures is real or an artifact of the underlying approximations (e.g., limited range of interaction parameters). The difficulties of converging a CE for multicomponent alloys has been also recently pointed out by Widom 16 on the example of finding the ground state of bcc NbMoTaW 17 where a conventional CE fails.
In the present work, we propose an alternative, highly efficient approach by fitting an accurate active-learning machine-trained potential which is used in subsequent Monte Carlo simulations.
Specifically, we employ a novel "on-lattice" machine-learning interaction model called low-rank potential (LRP) 18 , which is capable of including relaxation effects as well as it is capable of accurately representing interactions in multicomponent systems. It is thus well suited for dealing with a large number of components. As it was shown recently, the LRP requires fewer input structures to reach the same accuracy as a CE approach 18 . The model has only two adjustable parameters: an interaction cutoff radius and the approximation rank controlling the number of free parameters.
This allowed us to fit such models semi-automatically by active learning techniques. The approach can take local relaxation effects implicitly into account in a similar spirit as the CE (by allowing for local relaxations in the DFT calculations for the input structures, see Sec. 4). It also allows for a systematic estimation of errors in the predictions. We demonstrate the power of this new approach to efficiently and accurately explore huge configuration spaces by studying the finite-temperature phase stability of bcc NbMoTaW. Based on these studies we reveal a hitherto not reported chemical ordering at ambient temperatures. The impact atomic relaxations have on the phase stability and short-range order are discussed.

Results
An advantage of the approach based on machine-learning potentials, which will be introduced below, is that similarly to the CE formalism the configurations used to fit the potential can be chosen from static or fully relaxed (i.e., including local relaxation effects) calculations. In this way one can straightforwardly "switch on" and "off" the impact of local distortions and study their impact on phase stability and short-range order parameters.
In order to construct a machine-learning potential, we first generate a DFT dataset for fitting the potential as shown in the proposed workflow in Fig. 1. The initial training set consisted of 200 randomly generated configurations for which ground state DFT calculations were performed, each configuration constructed from a 2x2x2 bcc supercell (with 16 atoms) with randomly distributed atomic species. An initial low-rank potential (LRP) 18 is then fitted on the initial training set (details given in Sec. 4). In a nutshell, the LRP assumes a partitioning of the energy into contributions of each atomic environment, e.g., for a bcc lattice each environment is defined by nine atoms, i.e., a central atom and its eight nearest neighbors. This partitioning is different from those used in other formalism such as CE, which assumes a partitioning of the energy into individual two-body, three-body, etc., clusters.
For a given random 4-component alloy there are, in principle, 4 9 ≈ 250 000 possible atomic environments, however, a model with 250 000 free parameters is impractical to fit. In order to reduce the number of fitting parameters, the LRP makes a specific low-rank assumption on the representation of the interaction 19 (see Sec. 4). By applying this assumption we find that only about 500 parameters were eventually required to fit the potential with an accuracy of 1 meV/atom.
Although 250 000 environments are prohibitively too many for the purpose of fitting the energy contribution of each environment independently, we can still precompute and store the fitted data for each environment, which requires eventually only 2 megabytes of storage. This leads to an extremely fast energy evaluation which allows us to fit not only one, but an ensemble of 10 different potential models. An advantage of having such an ensemble of independent potentials is that not only interaction energies can be predicted, but the ensemble further enables us to estimate the predictive error (i.e., uncertainty) due to approximation of DFT energies with LRP in quantities of interest (such as formation energies) by observing the deviation of the different models. To be precise, we calculate the 95% confidence interval and call it the model uncertainty of LRP.
Based on the initial training set we evaluated the accuracy and performed LRP-based Monte Carlo simulations employing a larger 4x2x2 supercell including 32 atoms. From these simulations 100 additional, low-energy configurations were chosen, recalculated with DFT and added to the training set. These configurations were about 40 meV lower in energy than the randomly chosen initial training configurations, which corresponds to a thermal energy equivalent to ≈ 500K, a temperature close to the relevant phase transitions reported previously. This procedure has been repeated until the trained potential reaches a predictive error of less than 1 meV/atom. Our retraining (active learning) of the potential also ensures that the configurations span a wide range of compositions, from binaries to quaternaries as well as from chemically ordered to disordered configurations.
We first exemplify the performance of our proposed approach by training our potentials on a set of static calculations, i.e., with atoms dispersed on the ideal bcc positions not including the effect of local lattice relaxations. A set of 10 different potentials are trained and used in Monte Carlo simulations for a 12x12x12 supercell to compute the configurational contribution to the heat capacity, C V (T ), and the Warren-Cowley short-range order parameters 20 . Note that lattice vibrations not included in the present approach also contribute to the heat capacity. The results for C V are shown in Fig. 2 (blue solid line) and the predictive error is indicated in orange. First we note that the LRP model uncertainty is negligible for a large temperature range. Only close to the two observed phase transitions significant uncertainty due to the different LRPs is present. Note that we have chosen the number of MC steps sufficiently large so that the MC convergence errors are negligible, see Sec. 4. In order to evaluate our approach we first compare the results to previous works which were based on perturbational theories such as the generalized perturbation method 10 and also employed an ideal lattice approximation. In agreement with ref. 10 we find two phase transitions that the solid solution undergoes with decreasing temperature. The first transition when cooling down corresponds to a transition from the solid solution to a B2(Mo,W;Ta,Nb) ordered state as sketched in Fig. 2. We would like to emphasize that this B2-ordering has also been reported in all previous works so far 10-13, 15, 17, 21, 22 . Consistent with ref. 10 a second transition at low temperatures is identified as a phase decomposition into B2(Mo,Ta) and B32(Nb,W) (also sketched in Fig. 2).
The 30% decrease in the phase transition temperatures predicted by LRP as compared to 10 might also affect other properties such as the predicted degree of chemical short-range order at elevated temperatures. To elucidate this we focus next on the high-temperature short-range order where p ij is the probability of finding a j-type atom among neighbors of i-type atoms and c i , c j are the alloy concentrations of the corresponding elements. The results for the different α ij are shown in Fig. 3. Indeed, the supercell based calculations (dashed lines) reveal overall a smaller degree of short-range order at a given temperature as compared to the ones derived from the perturbation method (open diamonds). We find, however, that the main deviations are caused by the lower transition temperatures predicted by the presently used supercell based approach whereas the overall temperature dependencies of the α ij 's are rather similar. We note that the SRO for the Mo-Ta bonds is an order of magnitude larger than that for the other bonds as also observed in 10,15 . The question remains whether the assumption of static atomic positions (neglect of lattice relaxations) in these results causes any further qualitative change in the SRO predictions.
In order to account for local lattice distortions, we allowed for ionic relaxations in the configurations entering the training set. As we focus in the present work on the stability of the solid solution (homogeneous disordered alloy) we kept the cell volume and shape fixed to closer mimic the dominant randomized configurations. Note that if our objective was the zero-temperature ground state search, where phase decomposition and thus strong macroscopic volume fluctuations are anticipated, it might be important to consider volume relaxations. As mentioned above, the advantage of the employed LRPs is that it is straightforward to include local lattice distortions into the model.
We hence fitted a new ensemble of 10 potentials to the DFT data based on relaxed configurations and recomputed the SRO parameters.
The results for the newly trained LRPs including relaxation effects for the SRO are also shown in Fig. 3  Below 200K an ordered ground-state structure is formed consisting of 100 atomic planes consisting of the same type of atoms. The planes are repeated with the period of four lattice constants in the following sequence: Nb-Mo-Ta-W-W-Ta-Mo-Nb as illustrated in Fig. 4. We stress again that the calculations without relaxations resulted into a phase decomposition into B2(Mo,Ta)+B32(Nb,W) as also reported in ref. 10.
We next analyzed the MC-predicted structures below the first phase transition temperature.
We find that the solid solution turns not into a B2 ordered structure (as suggested in all previous works 10-13, 15, 17, 21, 22 and our non-relaxed scenario) but into a semi-ordered layered configuration as sketched in Fig. 4. This new state, which has not been reported before, can be characterized by atomic planes along the 100 direction occupied predominantly by one type of atoms, see At very low temperatures, a layered ordered structure, not reported so far, has been found. Interestingly, Widom has recently revisited bcc NbMoTaW and performed a very careful ground state analysis, and suggested a decomposition into a hR7(Mo 2 NbTa 2 W 2 ) and a cI2(Nb) structure 16 .
According to our calculations, its formation enthalpy is indeed by about 2 meV/atom below the low-energy structure we predicted if we allow, in addition to the ionic relaxations, for volume changes. This stabilization is mostly driven by pure Nb which has a much larger equilibrium volume than that of the solid solution. If we perform local relaxations only while keeping the volume fixed, the decomposition into hR7(Mo 2 NbTa 2 W 2 ) and cI2(Nb) is about 6.5 meV/atom above the here-found, layered ground state. This also highlights the performance of our approach in predicting the ground state under given external conditions, e.g., at a fixed lattice constant. We would like to emphasize that in contrast to the ground state considerations at zero K, the alloy remains macroscopically homogeneous under the phase transition occurring at room temperature.
We therefore do not expect that the inclusion of volume fluctuations would qualitatively alter our results at ambient temperatures.
In summary, we have proposed a new computational approach for the investigation of thermodynamic properties of high-entropy alloys. This approach is based on the LRP 18 , a computationally efficient machine-learning interatomic potential capable of accurately representing interactions in a system with many chemical components. The potentials are trained on DFT supercell calculations and thus allow to systematically include the impact of local lattice distortions. The approach is validated by employing a static-lattice approximation and comparing the results to existing ones 10,12,15 .
By including local atomic relaxations we found that, contrary to previous works, the solid solution is stable down to room temperature and transforms into a newly found, layered semi-ordered In the on-lattice model, LRP 18 , the energy of each configuration is partitioned into contributions of the separate atomic environments of each atom as where Ω is the lattice periodically repeated in space, V is the so-called interatomic potential defining the contribution of the atom at the lattice site ξ to the total energy E, σ(ξ + r i ) is the type of i-th neighbor of the atom located at ξ, r i is the vector connecting ξ with its neighbor, and n is the number of neighbors (including the central atom) that depends on the cutoff parameter. For bcc with only nearest neighbors included in the interaction, n = 9.
V can be thought of as a tensor with m n parameters, where m is the number of atomic types. In order to reduce the number of parameters, the low-rank tensor-train assumption is applied 19 . The tensor-train decomposition of V of rank r is simply where A i are matrices of rank r or less that depend on σ i ∈ {Mo, Nb, Ta, W}. The matrix A 1 has the size 1 × r, A 2 has the size r × r, etc., and A n has the size r × 1, so that the matrix product results into a scalar. 1 The matrix entries of A i are the parameters found from data. In this work we found that with r = 5 we can reach the accuracy of 1 meV/atom. This results in about 500 independent parameters.
There can only be 4 n ≈ 250 000 possible combination of input parameters (σ 1 , . . . , σ n ) of V . We therefore can precompute and store the values of V for all possible inputs. When this is done, the potential V is symmetrized over all 48 permutations of (σ 1 , . . . , σ n ) corresponding to the physical symmetries. This avoids artificial breaking of symmetry in the MC simulations and slightly improves the accuracy.
In order to find the parameters of V we minimize the following functional: where K is the number of the atomic configurations, σ (k) , in the training set (k = 1, ..., K), E Ä σ (k) ä is an energy of configuration predicted by the LRP model, E qm Ä σ (k) ä is the DFT reference energy.
The optimization functional Eq. 2 is not linear (because V depends nonlinearly on its tensortrain parameters). Therefore, there are plenty of local minima (by energy) in the space of the parameters. Thus, the minimization algorithm can find different local minima, depending on random initial parameters at the training stage. Therefore, thanks to the fast evaluation of an ensemble of 10 different LRPs were trained. Analysis of the independent predictions of the LRP ensemble makes it possible to estimate the uncertainty of the approach 18 .
The workflow of our calculations is shown in Fig. 1 For unbiased averaging, we employed the so-called burn-in 31 : we discard half of the MC iterations, and perform averaging over the second half. We note that discarding half of the iterations is necessary only near phase transitions, however, in order to keep our algorithm robust we always discard half of the iterations.
The most challenging MC simulations, for calculating the long-range pair correlation function at T = 240 K shown in Fig. 5, were conducted in a 48x48x48 cell. About 10 12 iterations were required to converge the graphs for the distance of up to 24 lattice constants, however, this was not enough to see the lack of long-range order perpendicular to the planes, in the [100] direction.
Calculations in even larger cubic cells, however, are not practical because of a very fast growth of the equilibration time with the cell size. Therefore, in order to see the lack of long-range order in the [100] direction we conducted calculations in a 100x20x20 cell. Note that these calculations do not provide any new information about the long-range order along the planes, but they reveal the long-range decorrelation perpendicular to the planes.