Short-range order in face-centered cubic VCoNi alloys

Concentrated solid solutions including the class of high entropy alloys (HEAs) have attracted enormous attention recently. Among these alloys a recently developed face-centered cubic (fcc) equiatomic VCoNi alloy revealed extraordinary high yield strength, exceeding previous high-strength fcc CrCoNi and FeCoNiCrMn alloys. Signiﬁcant lattice distortions had been reported in the VCoNi solid solution. There is, however, a lack of knowledge about potential short-range order (SRO) and its implications for most of these alloys. We performed ﬁrst-principles calculations and Monte Carlo simulations to compute the degree of SRO for fcc VCoNi, namely, by utilizing the coherent-potential approximation in combination with the generalized perturbation method as well as the supercell method in combination with recently developed machine-learned potentials. We analyze the chemical SRO parameters as well as the impact on other properties such as relaxation energies and lattice distortions.


I. INTRODUCTION
Multicomponent alloys, also known as high entropy alloys (HEAs) or chemically complex alloys (CCAs), have attracted enormous attention in the last decade, from theory and experiment, due to their remarkable materials properties and overwhelming compositional phase space for alloy design [1][2][3]. A recent example is an equiatomic face-centered cubic (fcc) single-phase VCoNi alloy revealing remarkable strength [4].
One of the key components in the design and exploration of multicomponent alloys is the knowledge about its phase stability. In fact, many such alloys have been demonstrated to decompose at elevated temperatures, e.g., the equiatomic fcc FeCoNiCrMn (also known as Cantor alloy) [5]. An important issue to address is therefore whether the solid solutions involve local chemical ordering or short-range order (SRO), which could affect, e.g., defect properties and hence their mechanical properties. This has, for example, been demonstrated in recent simulations of computed stacking-fault energies in fcc CrCoNi [6,7]. There are also some experimental results suggesting the presence of SRO in selected fcc HEAs such as, e.g., fcc CrCoNi [7,8], fcc FeCoNiCr [9,10], in Ni-Co-Fe-Cr * f.h.w.kormann@tudelft.nl alloys [11], and fcc CrFeCoNiPd [12]. In Ref. [7], a correlation between the mechanical behavior and SRO has been established for fcc CrCoNi.
For VCoNi, which is in the focus of the present work, a single-phase solid solution has been experimentally confirmed at 900 • C [4]. For lower temperatures, recent experiments revealed formation of a (Co,Ni)3V κ phase at around 850 • C [13] as well as a secondary, V-rich σ phase at temperatures below 800 • C. The κ phase features a close-packed ordered nine-layered hexagonal structure with an abcbcacab packing sequence [13][14][15][16] with characteristic layered-wise enrichment and depletion of V. This suggests that V could play a key role for the SRO in the solid solution. As mentioned above, the SRO could impact different properties such as, e.g., the degree of lattice distortions, which was one of the key computational descriptors to develop this high-strength alloy and for which V played a dominant role [4]. It is therefore crucial to accurately determine the degree of SRO and its implication for these properties.
As mentioned above, for most of the alloys only limited experimental information is available for quantifying the degree of SRO. In addition to experiments, in particular, ab initio-based density-functional theory (DFT) calculations have therefore been utilized to explore multicomponent alloys [2]. A number of such techniques are available to address phase stability and, in particular, SRO; among them the coherent-potential approximation (CPA) [17,18] based methods such as the generalized perturbation method (GPM) [19][20][21] in combination with Monte Carlo (MC) simulations or the concentration wave method [22,23]. These methods have been extensively employed in the past few years for multicomponent alloys [10,[23][24][25][26][27]] also due to their computational efficiency. A limitation of such mean-field approaches is, however, the great challenge to include local relaxation effects, which can be important for phase stability considerations [28,29]. One alternative is explicit DFT-based Monte Carlo simulations employing supercells [6,[30][31][32]. These calculations are, however, restricted to rather small simulation boxes (typically a few hundred) and also computationally extremely demanding. Therefore, DFT-derived energies are often mapped onto effective models such as the cluster-expansion technique [33], which is, however, computationally inefficient when treating multicomponent alloys [34]. An alternative to these approaches is provided by recently developed machine-learned potentials, which allow one to extremely efficiently parametrize the energy landscape, in particular for multicomponent alloys [29,34].
In the present work we employ two state-of-the art approaches. One is the recently developed low-rank potential (LRP), a machine-learned potential based on supercell calculations allowing us to fully account for lattice-relaxation effects. In addition, to estimate the importance of magnetic excitations, we perform Green's function multiple-scattering theory calculations to extract chemical effective interactions and study their dependence on the high-temperature magnetic state, which accounts for longitudinal spin fluctuations (LSF). The order-disorder transition is computed as well as the SRO parameters above the critical temperature. Potential implications, e.g., on lattice distortions and solid-solution strengthening contributions, are discussed.

A. Density-functional theory calculations for supercell approach
As outlined in the introduction above, the first model we employ is the low-rank interatomic potential (LRP) [34] (see also below), which is used as an interaction model in canonical Monte Carlo (MC) simulations in order to investigate short-range order and order-disorder phase transitions in the equiatomic VCoNi system. The LRP model is based on spin-polarized DFT calculations. To compute the reference energies of atomic configurations for the LRP training, VASP 5.4.1. [35][36][37][38] was used in combination with the projector augmented wave (PAW) method [39] and utilizing the Perdew-Burke-Ernzerhof generalized gradient approximation (PBE-GGA) [40]. Since we are interested in quantifying the degree of SRO at high temperatures, we performed the corresponding calculations at a lattice constant corresponding to high temperatures. The lattice constant has been derived from a Debye-Grüneisen model [41] at 1175 K, which is around the experimentally observed phase transition [13]. For the Debye-Grüneisen model, we computed the total energy of the alloy, E (V ), bulk modulus, and its derivative utilizing a 3 × 3 × 3 (108 atom) special quasirandom (SQS) [42] structure. The predicted room-temperature lattice constant of 3.599 Å is very close to the experimental value of 3.601 Å [4]. The obtained value at 1175 K, 3.643 Å, has been employed to compute the structures entering the training set, see below.
For the training set of the potential, different supercell sizes are chosen, ranging from 2 × 2 × 2 (32 atoms) to 3 × 3 × 3 (108 atoms). The value of the plane-wave cutoff energy is chosen to be 320 eV, which is 1.2 times larger than the recommended cutoff energy of the PAW pseudopotentials utilized. A k-point mesh is generated by the Monkhorst-Pack [43] scheme and is 6 × 6 × 6 for the system with 32 atoms, 4 × 6 × 6 for the system with 48 atoms, and 4 × 4 × 4 for systems with 108 atoms. Ionic relaxations of atomic positions are included in the calculations to account for the impact of lattice relaxations, which is presumably large in this alloy [4]. The energy convergence criteria for the ionic relaxations are set to achieve an error of 10 −4 eV. Electronic excitations are included by utilizing the Fermi distribution in the DFT calculations with a smearing parameter of 0.1 eV, which corresponds to a temperature close to the experimentally observed transition temperature.

B. The low-rank potential method
The atomistic structure for the LRP [34] model is represented as site occupancies in the ideal crystalline lattice. Each site is occupied by one of the atomic species (V, Co, Ni in this case). The sequence of atomic species among the closest neighbors of each site represents the environment of this site. The environment of each site has a contribution to the energy of the atomic configuration given by the formula (1) where V is the LRP model in the tensor form, as described below, ξ is the position of a central atom, σ (ξ + r i ) is the species of the ith site and r i is the vector connecting the central site with the ith neighbor, and n is the number of closest neighbors in the environment, including the central atom (n = 13 in the fcc case). The local lattice distortions are allowed as long as the topology of the supercell stays the same. The full energy of an atomic configuration can be described as a sum of individual contributions of the environments: where is the lattice sites periodically repeated in space. The tensor V contains energy contributions of all possible atomic environments, namely, m n environments, where m is the number of species in the system. Thus, for a fcc system with m = 3 species ({V, Co, Ni} in this case) the tensor V consists of about 1.6 million coefficients which is completely unfeasible to obtain from quantum-mechanical data. To remedy this, the LRP potential postulates a certain low-rank structure of V , or, to be precise, that there are some good approximants of V among the low-rank 13-dimensional tensors in the sense of the tensor-train decomposition [44]. This decomposition scheme will be explained in the following. Let us consider for the moment that V is a large N × N matrix (i.e., a twodimensional tensor). We identify the matrix V with a discrete function V (σ 1 , σ 2 ). Assuming that the matrix V has a low rank r means that the column vectors of V can be expressed as a linear combination of r vectors w 1 , . . . , w r : Here we used the notation u i (σ 1 ) for the coefficients of the linear expansion to make the notation symmetric. The 113802-2 expansion (3) is reminiscent of the method of separation of variables and could be motivated as follows: Often when the matrix V describes some regular (smooth) function, only a few singular values s i are far from zero, where the right-hand side is precisely a singular value decomposition. For symmetric positive-definite matrices V , s i would be the nonzero eigenvalues with u i = w i being the corresponding eigenvectors. Adsorbing s i into u i and w i yields (3). One could generalize (3) to more than two dimensions by postulating V (σ 1 , . . . , σ n ) = r i=1 a 1 (σ 1 ) · · · a n (σ n ). This formalism leads, however, to computational issues when used in optimization algorithms [45]. Instead, we utilize the densitymatrix renormalization-group analogy, and hence cast (3) in the vector product form: with a 1 and a 2 being some r-dimensional vectors. We then generalize: where A i are matrices with rank r or less (A 1 , A n are of size 1 × r and r × 1, correspondingly, A 2 , . . . , A n−1 are of size r × r), r is this rank of the decomposition. The elements of the matrices A i are the parameters of the LRP. Thus, the tensortrain decomposition reduces the number of these parameters from m n to nmr 2 . In the context of our three-component alloy, σ i can take either of the three values (V, Co, or Ni). We found that r = 5 is an optimal rank and thus about 900 parameters were fit instead of 1.6 million. The parameters are found as a result of solving the minimization problem with the following functional: where σ (k) are the atomic configurations, the total number of which in the training set is K, and E (σ (k) ) and E qm (σ (k) ) are the energies of σ (k) predicted by LRP and calculated by DFT, correspondingly. The minimization was done with an alternating least squares method that reduces to simply optimizing one matrix A i at a time (this is a linear problem), and simulated annealing which consists of adding random Gaussian noise to every element of A i with a magnitude that gradually reduces from one iteration to the next. The dependence of V on its parameters is not linear. Thus, different local energy minima exist in the parameter space. Therefore, depending on the initial parameters, the minimization algorithm finds different local minima. The latter means that each LRP gives independent energy predictions, and with a trained ensemble of several LRPs, the uncertainty level of the LRP model can be estimated. In this work, an ensemble of 10 LRPs was used. The 10 LRPs differ by the random initialization of the matrices A i in Eq. (4). The coefficients of each matrix were drawn from a normal distribution with zero mean. The value of the variance does not affect the results because, at each step of the optimization, the matrices A i are rescaled to the data [34]. We find that, in practice, all 10 LRPs have comparable fitting errors.
The initial training set consisted of random disordered configurations with sizes 3 × 3 × 3 (108 atoms), 3 × 2 × 2 (48 atoms), 2 × 2 × 2 (32 atoms) confined to equimolar (for 108-and 48-atom cells) or as close to equimolar as possible (for the 32-atom cell). The number of configurations with 108, 48, and 32 atoms was 86, 70, and 70, correspondingly, in the training set, and 10, 10, and 9 configurations in the validation set. We note that fcc point symmetries were applied, so that each configuration enters the training set together with the other 47 equivalent configurations, thus effectively generating about 10 000 data points to determine the 900 parameters of the LRP.
To improve the accuracy of the LRP, the workflow from Ref. [29] (see Fig. 1) for selecting new configurations is used. MC calculations are conducted for systems with 108 atoms. The values of the enthalpy and heat capacity for the LRP ensemble are compared for different temperatures, and configurations that correspond to the temperature range with the highest deviations are sampled.
The sampled configurations are added to the training and validation set and the new ensemble of LRP is trained. Each training starts with a new random distribution of initial matrices A i . This procedure continues until the training error stops changing significantly. In this work, 40 new configurations were sampled and added to the data sets.
The rank is an adjustable parameter and, to choose it, we considered values from three to six. Starting from r = 5, the accuracy of prediction achieves 2 meV/atom.
LRP is a representative of the machine-learned class of interatomic potentials and can be compared with the cluster expansion method (CE). In contrast with the CE [46], LRP has a form of many-body interaction without separating it into two-body, three-body, etc. Similarly to CE, if atomic displacements do not change the topology of the lattice structure, then the local lattice distortions can be taken into account.

C. Monte Carlo method
The LRPs are used in a canonical Mote Carlo method [47]. We mainly focus on the temperature range of 1100-1500 K, 113802-3 as the experimentally observed order-disorder phase transition is at approximately 1175 K [4].
The simulations are carried out for systems with 6912 atoms (12 × 12 × 12 in lattice units, based on a four-atom primitive fcc cell). The number of MC steps is adaptive for different temperature ranges: for temperatures higher than 2000 K-2 × 10 8 steps, for the interval of 1250-2000 K-2 × 10 9 steps, for temperatures lower than 1250 K, the number of steps is 2 × 10 10 . For the smaller structures with 32 atoms (2 × 2 × 2 in lattice units), the number of steps is two orders of magnitude smaller. To achieve an unbiased averaging, the so-called burn-in approach [48] was implemented, i.e., we neglected the first half of MC steps for each temperature value. This technique is necessary at the range of temperatures that include phase transitions but was used permanently in our calculations to improve the robustness of the algorithm.

D. Multiple-scattering theory calculations
As a second approach, in order to obtain a qualitative picture of the impact of low-and high-temperature magnetism as well as the effective interactions in V-Co-Ni alloys, the Green's function multiple-scattering theory was used as implemented in the exact muffin-tin orbital method (EMTO) [49][50][51]. The electronic structure of the random alloys was obtained within the coherent-potential approximation (CPA) as well as by using the locally self-consistent Green's function (LSGF) technique [52,53] within the EMTO method, ELSGF [54].
All the EMTO-CPA calculations were done by utilizing the Lyngby version of the Green's function EMTO code, which allows us to calculate effective interactions [55] employing the screened generalized perturbation method (SGPM) [19,56,57], to include the impact of longitudinal spin fluctuations (LSFs). Moreover, the screened Coulomb interactions are treated properly in the single-site DFT-CPA approximation [20,21].
The self-consistent electronic structure calculations were performed in the local density approximation (LDA) using the Perdew and Wang functional [58] while the total energies were calculated by the full charge-density technique [51] employing PBE-GGA [40]. The Brillouin-zone integration was done by using a 30 × 30 × 30 Monkhorst-Pack grid [43] or finer. All calculations have been done with l max = 3 for partial waves, and the electronic core states were recalculated at every iteration during the self-consistent calculations for the valence electrons.
The ELSGF method was used to calculate the screened Coulomb interactions used in the DFT-CPA part of the EMTO-CPA calculations. In this case, the one-electron potential of the alloy components and the total energy have additional contributions, V i scr , and E scr , respectively, due to the screening charge around atomic spheres, which is not accounted for in the single-site approximation [20,21]: Here,q i and α 0 i are the net charge of the atomic sphere of the ith component in the single-site CPA calculations and its onsite screening constant, respectively; S WS is the Wigner-Seitz radius; E i scr is the contribution of the screened Coulomb interactions to the electrostatic energy of the ith alloy component; β scr is the average on-site screening constant, which accounts for the electrostatic multipole moment energy contribution due to the inhomogeneous local environment of different sites in a random alloy.
The screening constants for the Co-Ni-V alloys were determined from ELSGF 864-atom 6 × 6 × 6 supercell calculations. The on-site screening constants α 0 i were determined from the conditional average of the net charges q i and the Madelung potentials V Mad The intersite screening constants α i j p are needed in the calculations of the electrostatic contribution V i j−scr p to the SGPM potential at the pth coordination shell for the i-j pair of alloy components, Here,q i j =q i −q j were obtained in the supercell ELSGF calculations for random alloys from the screening charge by exchanging the corresponding alloy components (i and j in some specific sites have random local environment on average) as described in Refs. [20,21].

A. Phase stability and short-range order
A goal of the present work is to investigate the degree of short-range order in VCoNi alloy as well as the interplay with lattice distortions and relaxation effects, which is closely related to the exceptional strengthening of the VCoNi alloy [4].
As sketched in Fig. 1, we first included 226 configurations to the training set. We trained an ensemble of 10 LRPs, which were used in the MC simulations for systems with 108 atoms. We analyzed the dependencies of enthalpy and heat capacity on temperature obtained with these 10 LRPs. Based on this initial set of potentials, we sampled a number of new configurations at a range of temperatures close to the precomputed order-disorder phase transition, i.e., near temperatures where significant fluctuations among the initial potentials was observed. As a result, 40 new configurations were selected and included in the fitting procedure, namely, 10 configurations for each of the following temperatures: 760, 860, 1100, and 1140 K. These configurations were recalculated with VASP (see above) and added to the training and validation sets, thus improving the prediction error of the finally utilized LRPs down to 2 meV/atom.
We also evaluated the impact of including a number of DFT-computed 4 × 4 × 4 supercell results to the training set. The accuracy of the energy prediction is still about 2 meV/atom, resulting in qualitatively similar predictions of, e.g., specific-heat capacities, as will be shown below. We thus concluded that the inclusion of explicit DFT supercell calculations beyond the mainly used 3 × 3 × 3 cells into the LRP training set would not qualitatively alter our main results and conclusions.
To investigate the temperature-dependent evolution of SRO and phase stability of the alloy, LRP-based MC simulations were carried out for larger 12 × 12 × 12 supercells (6912 atoms). In Fig. 2 the dependence of specific-heat capacity on temperature, C V (T ), is presented. We first observe a phase transition with a characteristic peak in C V (T ) around 1400-1500 K. There is a slight dependence of the transition temperature depending on the largest DFT-computed supercells in the training set, which is manifested in deviations mainly around the peak whereas the lower and higher temperatures appear less size-dependent. The overall transition temperature is also somewhat higher than the experimental one, which is between 1123 and 1173 K [4,13]. The discrepancy may be related to missing thermal fluctuations such as, e.g., vibrations which could in principle alter the predictions [59]. We also note that the present approach cannot account for the experimentally observed hexagonal κ phase. On the fcc lattice we observe a phase transition to a partially ordered structure, where V atoms almost fully occupy one of the four primitive cubic sublattices of the fcc structure, as shown in the inset of Fig. 2. The computed atomic concentration of V for each of the four sublattices are shown in Fig. 3. This clearly shows that the observed phase transition is caused by a strong site preference of V on one of the four sublattices, while the remaining V is equally distributed on the other sublattices, which corresponds to the M3V-L1 2 ordering.
The L1 2 structure is formed in Co 3 V alloys as a metastable phase if the alloys are quenched and constrained to the fcc lattice (see, e.g., Ref. [60] and references therein). Recently, a similar L1 2 ordering has been also reported for a fourcomponent fcc FeCoNiV alloy [61]. It may therefore be that the L1 2 ordering observed in the present work for VCoNi can exist as a metastable phase. The experimentally observed σ phase [4,13,62] at temperatures below 800 K cannot be covered by our simulations because we are operating on the fcc lattice. We therefore focus in the following on the SRO effects at high temperatures.
To quantify the degree of SRO, the Warren-Cowley SRO parameters [63] were calculated: where α i j m is the Warren-Cowley SRO parameter for the atomic types i and j at the mth coordination shell; p i j m is the probability of finding atom type j at the mth coordination shell i of atom m, and c i , c j are the concentrations of elements i and j in the alloy. According to this definition, positive (negative) values of the SRO parameter at the mth coordination shell mean that atoms i and j avoid (attract) each other at the corresponding coordination shell.
The computed SRO parameters for the first and second coordination shells are shown in Figs. 4(a) and 4(b). Obviously, there is strong ordering between pairs Co-V and Ni-V at the first coordination shell, while these pairs exhibit strong repulsion at the second coordination shell, which corresponds to the L1 2 type of ordering (in the completely ordered A 3 B-L1 2 phase, the Warren-Cowley SRO parameters at the first two coordination shells are −1/3 and 1). Thus, the formation of atomic SRO in VCoNi alloys is mainly driven by the strong interaction of V with Co and Ni.

B. Interplay and impact of short-range order and lattice distortions
Since the fcc VCoNi solid solution is prone to significant local lattice distortions [4], it is important to evaluate how relevant the inclusion of the corresponding relaxation energies are for an accurate modeling. An advantage of the LRP approach employed here is that local distortions can effectively be "switched off," allowing us to explore their impact and relevance, e.g., for the prediction of SRO parameters as well as phase stability.
Similarly to our previous work [29], we evaluate the impact of local distortions and the corresponding relaxation energies by training a new ensemble of LRP potentials on a set of static DFT calculations followed by subsequent MC runs, as discussed above. Thus-obtained specific-heat capacity is shown in Fig. 5 and shows that the M3V-type of ordering sets in only above temperatures of 2000 K. This is about 30% larger compared with the case when relaxation effects are included and highlights the importance of taking relaxations into account for modeling such alloys.
To further elucidate the impact of relaxation energies, we performed an additional MC simulation from the previous LRPs, including relaxations for a smaller 4 × 4 × 4 supercell (containing 256 atoms) and selected ten snapshots at three representative temperatures: at 1000 K for a M3V type of ordered structure, at 1540 K for a strongly SRO-containing solid solution, and at 4000 K to represent a random solid solution. In total, 30 new structures were recalculated with DFT to extract the relaxation energies, which are depicted in Fig. 6 (light blue bars). The error bars indicate the standard deviation for the mean value as obtained from averaging over the ten individual calculations for each case. There is a clear trend of increasing relaxation energies with increasing disorder. The largest relaxation energies of roughly 30 meV/atom are found for the solid solution and are almost twice as large as found for the M3V type of ordered structures. This is consistent with our findings that the computed transition temperature decreases drastically if relaxation effects are included in the simulations.
Based on the examples considered, we also evaluated the impact of the electronic free-energy contributions to the phase stability. As mentioned above, all our calculations include electronic free energies based on a Fermi smearing parameter of 0.1 eV, which corresponds to a temperature of 1160 K, close to the experimentally observed phase transition. In Fig. 6, the electronic free-energy contributions are shown in light green for the different scenarios. Although the electronic free energies are overall of similar magnitude as compared with FIG. 6. Electronic free-energy contributions at 1160 K (corresponding to a Fermi smearing of 0.1 eV) in light green, relaxation energies (light blue), and mean-squared atomic displacements (light red) obtained from DFT calculations of 256-atom (4 × 4 × 4) supercells selected from MC runs for three different scenarios: MC simulations at 1000 K (M3V type of ordering), at 1540 K (SROcontaining solid solution), and at 4000 K (random solid solution). Relaxation energies and lattice distortions increase with increasing disorder. The standard deviation obtained by averaging over 10 different supercells for each case is depicted as error bars. the relaxation energies, the relative changes are much smaller, implying that SRO does not significantly alter the electronic density of states, and vice versa, the computed phase transition is not largely affected by electronic excitations. There is overall a slight increase in the electronic free energy with increasing disorder, indicating that the electronic excitations contribute to the stabilization of the solid solution.
We next considered the reverse impact of relaxations and SRO, i.e., how largely SRO can impact the amount of local lattice distortions. For this purpose we resort to the mean square atomic displacement (MSAD), defined as where i runs over all atomic positions R i , R ideal i denotes the ideal fcc lattice positions, and N is the number of atoms i.
In previous works [4,64], the square root of MSAD had been utilized as an effective parameter to quantify the degree of lattice distortions as well as a promising descriptor correlating with the mechanical strength of the alloy. It is therefore important to evaluate how the distortions are affected by the degree of SRO found in the alloy. We therefore resort to the 30 DFT calculations at the three representative temperatures discussed above and evaluated the MSAD parameter for all of them. The results are also shown in Fig. 6 (light red colored bars).
Similar to the relaxation energies, we again find an increase in lattice distortions with increasing disorder. The relative change is, however, less than with the relaxation energies. For example, for the M3V type of ordered structures, the square root of MSAD drops by about 25%. For the strongest SROcontaining solid solution case chosen at 1540 K, just above the transition temperature (see Fig. 2), the decrease in the lattice distortion parameter is only 6%. Assuming that the MSAD value correlates with the alloys' solid-solution strengthening ability [4,64], it would suggest that SRO has only very little impact on this mechanism. Note that other defect energetics, such as stacking-fault energies, might be more sensitive to SRO as, e.g., recently discussed for CrCoNi [6,7].

C. Finite temperature magnetism and effective interactions in VCoNi alloys
The low temperature magnetism of VCoNi alloys depends on the alloy content and the state of order. Although fcc Co has the highest Curie temperature among all the transition metals, it is a relatively weak itinerant magnet, i.e., the local magnetic moment of Co atoms is quite sensitive to the magnetic state, the local structure, and the chemical environment of Co atoms. At 0 K, e.g., fcc random Co-V alloys are nonmagnetic if the concentration of V is larger than 50 at.%.
The ferromagnetic state of pure Ni is much weaker than that of Co, and Ni loses its magnetic moment in fcc random Ni-V alloys when the concentration of V exceeds 12 at.%. One would therefore expect the ordered Ni 3 V to be nonmagnetic. This is indeed true for the ground-state DO 22 structure. However, the L1 2 -Ni 3 V is ferromagnetic with the main contribution to the total magnetic moment of the alloy coming from the V atoms, which is on the order of 1μ B , while the magnetic moment of Ni is only 0.05μ B and is parallel to that of V.
Exchanging some Ni by Co (randomly) in L1 2 -Ni 3 V changes the ground-state magnetic structure: in the L1 2 -(Co 50 Ni 50 ) 3 V alloy, the magnetic moment of Ni remains small, of the order of 0.05μ B , and is parallel to that of Co, which is about 0.5μ B , while the magnetic moment of V is about 0.4μ B and it is antiparallel to those of Co and Ni. If one now increases the concentration of V towards the equimolar VCoNi composition in the L1 2 structure by placing randomly V atoms on the Ni and Co sublattices, the alloy becomes nonmagnetic. At the same time, the equimolar fcc VCoNi random alloy is ferromagnetic, with the following magnetic moments of Co, Ni, and V: 0.6μ B , 0.07μ B , and −0.1μ B , respectively. However, its magnetic energy, i.e., the difference of the total energy of this alloy in the ferromagnetic and nonmagnetic states, is only about 8 K/at. This is a short description of what happens with the local and total magnetic moments and magnetic state in V-Co-Ni alloys at 0 K, which have been obtained in the EMTO-CPA calculations in this work. It is clear that, independently of the ordered state, VCoNi alloys can be considered as nonmagnetic at very low temperatures. However, this cannot a priori be assumed at high temperatures, where entropy induces large longitudinal spin fluctuations in all three alloy components. In particular, the disordered local moment (DLM) calculations [65] with the LSF entropy 3 ln( m ) [66] show that, at 1000 K, the magnetic moments of Co, Ni, and V are 1.33μ B , 0.65μ B , and 0.87μ B , respectively. These are quite sizable magnetic moments stabilized by the LSFs and one can therefore expect that they may also affect the effective interactions in this alloy, which we discuss below.
Following Ref. [55], we consider here for clarity a quasibinary form of the pair interaction contribution to the configurational Hamiltonian: where V (2)-αβ p are the effective pair interactions of the α and β alloy components at the pth coordination shell; δc α i = c α i − c α is the concentration fluctuation of the α component from its average concentration c α i in the alloy at site i. In Fig. 7, we show the effective pair interactions for all the three pairs of alloy components in equimolar VCoNi alloy obtained from the SGPM nonmagnetic and paramagnetic DLM-LSF calculations. The DLM-LSF calculations were performed at 1000 K. These calculations are done for the ideal fcc lattice positions, without consideration of local lattice relaxations, which are quite significant in this system due to relatively large size mismatch of Co-V and Ni-V alloy pairs, as has been demonstrated above.
To get the total effective interactions, the strain-induced interactions of Co-V and Ni-V pairs has therefore been estimated in the dilute limit of pure Co and Ni considering V as an impurity. The strain-induced interactions of the Co-Ni pairs were neglected since the size difference of these alloy components is insignificant, at least compared with those of Co-V and Ni-V alloys [67]. These calculations were also done with VASP using 256-atom (4 × 4 × 4) supercells in the ferromagnetic state to mimic the increased size of alloy components at high temperatures due to the quite substantial local magnetic moments. Note that the strain-induced interactions are mainly due to the size effect, while SGPM catches the "chemical" part of interactions, related to the renormalization of the electronic states.
Obviously, the main effect due to strain-induced interactions and LSF is related to a renormalization of the effective interactions at the first coordination shell. The LSF effect is more pronounced in the case of Co-V interactions due to relatively large local magnetic moments of Co and V, while the local lattice relaxation effects are stronger in the case of Ni-V pairs, as Ni atoms are somewhat smaller than Co ones. The important point here, however, is that all these renormalizations only change the interactions quantitatively, not qualitatively. Both, Co-V and Ni-V interactions still promote the so-called (100)-type ordering specific for the L1 2 structure on the fcc lattice.
There is, however, an interesting observation if one compares the pair interactions obtained in the VCoNi alloy with those obtained for the binaries (shown in the insets in Fig. 7). In fact, the Co-V and Ni-V interactions in binary Co 3 V and Ni 3 V alloys are very similar (see red and black diamonds in insets in Fig. 7). The ground state of Ni 3 V is, however, DO 22 and not L1 2 . This suggests that the DO 22 structure is stabilized by multisite interactions, which is indeed the case. The multisite interactions are also quite strong in the equimolar VCoNi alloy, especially the three-site interactions. However, their effect on ordering in VCoNi is very small due to the equimolar, i.e., also more "symmetric," alloy composition because the three-site interactions contribute to the ordering only for asymmetric alloy compositions.
Thus, the pair effective interactions play the most important role for atomic ordering in the equimolar VCoNi alloy, and the Monte Carlo results presented in Fig. 4 can be readily understood from the pair effective interactions in Fig. 7: strong ordering of Co-V and Ni-V atoms at the first coordination shell and at the same time a quite strong repulsion at the second coordination shell. The Co-V and Ni-V SRO parameters are extremely similar in the Monte Carlo simulations presented above, which is due to the fact that the effective interactions for these pairs of alloy components are very close to each other. And, finally, the atomic SRO parameters of NiCo are also consistent with the trend in the corresponding Ni-Co effective interactions in Fig. 7.

IV. CONCLUSIONS
The phase stability and short-range order (SRO) in the medium entropy fcc VCoNi alloy have been investigated by utilizing a recently developed machine-learned potential in combination with DFT supercell calculations as well as employing a complementary DFT-based multiple-scattering theory. On the fcc lattice, a phase transition to an M3V-type of ordering is found at about 1500 K. Above the phase transition, the SRO is mainly caused by a strong ordering of Co-V and Ni-V pairs at the first coordination shell accompanied by their relatively strong repulsion at the second coordination shell. By using two complementary DFT approaches, we showed that the impact of longitudinal spin fluctuations on the effective pair interactions is small compared with the chemical contribution and does not change the qualitative trends.
The inclusion of atomic relaxations and relaxation energies was found to be important for a quantitative description, and neglecting it results in an increase of the predicted transition temperatures by more than 30%. This is caused by the large relaxation energies of the VCoNi solid solution. On the other hand, the relative lattice distortions remain rather strong, also in the SRO-containing alloys. A small decrease in the mean square atomic displacement value suggests that SRO would affect the overall lattice distortions only very little, and hence likely the friction stress and solid solution ability. Finite-temperature electronic free-energy contributions also contribute to the stabilization of the random alloy, although their relative impact is weaker as compared with the effect of lattice relaxations.
(Linköping) and PDC (Stockholm) resources provided by the Swedish National Infrastructure for Computing (SNIC). A.V.R. also gratefully acknowledges the financial support under the scope of the COMET program within the K2 Center "Integrated Computational Material, Process and Product Engineering (IC-MPPE)" (Project No COMET 859480). This program is supported by the Austrian Federal Ministries for Climate Action, Environment, Energy, Mobility, Innovation and Technology (BMK) and for Digital and Economic Affairs (BMDW), represented by the Austrian research funding association (FFG), and the federal states of Styria, Upper Austria, and Tyrol.