Computational Materials Design of High-Entropy Alloys Based on Full Potential Korringa-Kohn-Rostoker Coherent Potential Approximation and Machine Learning Techniques

Kazunori Sato; Genta Hayashi; Kazuma Ogushi; Shuichi Okabe; Katsuhiro Suzuki; Tomoyuki Terai; Tetsuya Fukushima

doi:10.2320/matertrans.MT-MG2022012

Abstract

Computational materials design (CMD) based on the first-principles electronic structure calculations is demonstrated for two topics related to the design of high-entropy alloys (HEAs). The first one is a construction of prediction model of elastic constants. By applying machine learning technique with the use of the linearly independent descriptor generation method to the database of elastic constants of 2555 BCC HEAs generated by the full potential Korringa-Kohn-Rostoker coherent potential approximation (FPKKR-CPA) method. The obtained model is used to predict new HEAs with high Young’s modulus. The second topic is a simulation of atomic arrangement in HEAs at finite temperature. In this simulation, HEAs are described by using the Potts-like model and the interaction parameters are determined based on the generalized perturbation method combined with the KKR-CPA method. Monte Carlo simulations for the models of CrMnFeCoNi and CrMnFeCoCu predict atomic arrangements which are consistent to the experimental observations.

1. Introduction

High-entropy alloys (HEAs) are composed of more than 5 different elements with nearly equal atomic concentration.¹⁾ HEAs attract much attention due to their possibility to realize various functionality depending on the choice of the constituent elements. Previously reported functionalities include high strength at elevated temperature,²⁾ favorable temperature dependence of strength and ductility,³⁾ thermal stability,⁴⁾ corrosion resistance,⁵⁾ fatigue resistance,⁶⁾ bio-compatibility⁷⁾ and so on. Generally, the constituent elements are chosen from the following point of view, namely the high-entropy effect, lattice distortion effect, sluggish diffusion effect and cocktail effect.¹⁾ The effectiveness of these core effects is still under investigation and reliable guiding principles to propose new HEAs with desired functionality have not been established yet.

The first-principles electronic structure calculations might be a good choice to realize effective and practical method for designing functional HEAs, since a number of combinations of constituent elements is too large to allow exhaustive experiments. This is why many computational researches based on the electronic structure calculations by the density functional theory (DFT) have been reported so far.⁸^–¹¹⁾ However, due to the structural complexity and combinatorial explosion, which are inherent attributes of HEAs, straightforward application of the first-principles calculations to the computational materials design (CMD) of HEAs is still impractical.

As a practical solution to realize CMD for HEAs, recently an approach of materials informatics (MI) is focused. For example, the first-principles calculations are performed for a limited number of HEAs to generate a database and by employing machine learning techniques to the database one can obtain prediction models for physical properties of interest.¹²^,¹³⁾ In this paper, as one of such attempts, we demonstrate that the combination of the full potential Korringa-Kohn-Rostoker coherent potential approximation (FPKKR-CPA)¹⁴⁾ and the linearly independent descriptor generation (LIDG) method¹²^,¹³^,¹⁵⁾ works reasonably and effectively for the prediction of elastic properties of HEAs.

The above MI approach by the combination of FPKKR-CPA and LIDG can be a reasonable starting point of CMD for HEAs, but it is not enough to describe all of the complexity in real HEAs. For example, in the CPA calculations completely disordered configuration is assumed and the fluctuation of atomic arrangement such as phase separation or short-range order¹⁶^–¹⁹⁾ cannot be included. In this paper, by using the generalized perturbation method (GPM)²⁰⁾ we calculate effective pair interactions (EPIs)²¹^,²²⁾ between atoms and try to simulate atomic arrangement in HEAs at finite temperature by using the Monte Carlo simulation (MCS) to assess how realistic the assumption of complete disorder is in HEAs.

2. Methods

In this paper, we will show 2 different examples of computational materials design. The first one deals with a model construction for the materials properties of HEAs based on the first-principles calculations. As typical physical property, we focus on the elastic constants of BCC HEAs. The second one is a simulation of atomic configuration of FCC HEAs at finite temperature. Theoretical frameworks for these computations will be given in the following subsections.

2.1 First-principles calculations and model construction for elastic constants of HEAs¹³⁾

In order to construct a prediction model of elastic constants of HEA, firstly we have performed total energy calculations within the framework of the DFT by deforming cubic unit cell of HEAs systematically. We employ the FPKKR-CPA method¹⁴⁾ for the DFT calculations for the following two reasons. Firstly, by using the CPA method one can calculate configuration average of the electronic structure of disordered HEAs without using a large super-cell.²³⁾ Thus, the computational effort is greatly reduced. Moreover, after taking configuration average the cubic symmetry of the system is recovered and we only need 3 different deformations to determine 3 independent elastic constants. In the super-cell calculations, complicated averaging procedures might be required to obtain elastic constants of HEAs since each super-cell does not have cubic symmetry due to the random atomic configuration. Details of the FPKKR-CPA calculations for elastic constants of HEAs were already explained in our previous paper.¹³⁾ We simply use the database, which consists of bulk modulus B, two cubic share moduli c′, c₄₄ and Young’s modulus of 2555 equiatomic quinary BCC HEAs, generated in Ref. 13).

Concerning to the linear regression of calculated elastic constants, the descriptors are generated based on the following physical quantities of constituent elements such as 3 independent elastic constants (B, c′, c₄₄), lattice constant, group and period of the elements, atomic number and electron density parameter (r_s) of components of HEA. Here, r_s is defined as $\frac{4\pi }{3}r_{s}^{3} = \frac{1}{\rho }$, where ρ is the electron density in the interstitial region. From these basic quantities we generate linearly independent descriptors systematically. For the training process, 80% of the data are used and the remaining 20% is used only to test the obtained model. See Ref. 13) for a detailed description of descriptors and optimization procedure of the model.

2.2 Monte-Carlo simulation of atomic configuration in HEAs

For the assessment of the stability of disordered phase of HEAs at finite temperature, the HEAs is described by using the 5-state Potts-like model.²⁴⁾ The Hamiltonian is written down as, $H = - \sum\nolimits_{\langle i,j\rangle }^{N}V (\alpha ,\beta )_{ij}\delta (\sigma_{i},\alpha )\delta (\sigma_{i},\beta )$, where V(α, β)_ij is the EPI between site i and j. σ_i indicates the elements at site i and takes any integer from 1 to 5 which distinguishes the element occupying site i. δ is the Kronecker delta, thus corresponding EPI is chosen via α and β and sum up to get configuration energy. For element α and element β, the EPI between site i and j is defined as V(α, β)_ij = E(α, α)_ij + E(β, β)_ij − 2E(α, β)_ij, where E(α, β)_ij corresponds to the total energy when element α is located at site i and element β is located at site j. By definition, positive (negative) EPI indicates repulsive (attractive) interaction. The EPIs between constituent elements are calculated by the generalized perturbation method.²⁰⁾ For the present calculations we employed formulation given in Ref. 21), and EPI’s are calculated as $V(\alpha ,\beta )_{ij} = - \frac{2}{\pi }\textit{Im}\int^{E_{F}}dE\ Tr[\Delta_{i}(E)\tau_{ij}(E)\Delta_{j}(E)\tau_{ji}(E)]$, where Δ_i(E) = X_α(E) − X_β(E) and τ_ij(E) = ∑_k[t⁻¹(E) − g(k, E)]⁻¹ exp{ik·(R_i − R_j)}. X_α is a single site scattering matrix due to element α, t is coherent t-matrix and g is the KKR structure constant. τ_ij is called scattering path operator. E and k are energy and wave vector of electron, respectively. R_i is the lattice vector of site i. Detailed explanation and derivation can be found in Refs. 21) and 22). For the calculation of EPIs we employed the subroutine developed by Long et al.²⁵⁾ implemented to Akaikkr package.²³^,²⁶⁾ Once the EPI’s are determined we can generate atomic configurations in HEAs at finite temperature by performing MCS.

In performing the CPA, single site approximation is usually employed for calculating impurity Green’s function. This approximation causes calculation error in averaging Madelung energy especially when charge transfer between constituent elements is large. Fortunately, this error can be partly repaired by assuming short-range screening potential.²⁷^,²⁸⁾ For the present calculations, EPI’s were calculated by including such screening. Another error might come from local lattice relaxation. This error is inherent in CPA and the calculated EPI’s might have considerable error especially in cases where atomic sizes of constituent elements are so different.

For the MCSs we use homemade program code and simple Metropolis algorithm is employed. The MCS is performed for 10 × 10 × 10 super cell which is constructed from the conventional BCC unit cell. The initial configuration is prepared by assuming random distribution of constituent atoms in the super cell, and for each temperature the MC sampling was done 5000 times for every 5 MC steps for taking thermal average. The simulation starts from highest temperature (1,500 K), and calculate temperature dependence of the Warren-Cowley short range order parameter (SRO) with decreasing temperature by 100 K step. The present approach for the simulation of atomic configuration in disordered systems were previously applied to the semiconductor spintronics materials and succeeded to explain nano-structure formations in dilute magnetic semiconductors.²⁹^–³²⁾

3. Results and Discussions

3.1 Exhaustive search of functional HEAs

According to the method briefly described in the sec. 2.1, we have constructed the prediction model of B, c′, c₄₄ of BCC HEAs. From these 3 elastic constants, Young’s modulus is calculated. The obtained accuracy of the model is 10.2 GPa for B, 4.5 GPa for c′ and 2.4 GPa for c₄₄. For these models of B, c′ and c₄₄, 160, 208 and 198 descriptors are optimally chosen from 2nd order products of basic descriptors, respectively.¹³⁾

Based on the constructed prediction model, we tried exhaustive search of HEAs. Table 1 shows the list of high Young’s modulus HEAs predicted by the model. Our model does not predict the stability of BCC HEAs, thus we cannot judge the synthesizability of the predicted HEAs. In order to make the prediction as realistic as possible, we consider empirical rule¹^,³¹^–³³⁾ for the existence of HEAs when we list up the HEAs in Table 1. The detailed description of the criteria, based on topological parameter (δ), enthalpy of mixing (ΔH), thermodynamic parameter related to melting temperature (Ω), and valence electron concentration (VEC), can be found in the Refs. 1), 13), 33)–35). It can be argued that presently constructed prediction model for elastic constants of BCC HEAs might be useful for exploring large material space of HEAs.

Table 1 Prediction of HEAs with high Young’s modulus E by using the model constructed based on the first-principles calculations. In addition to E some parameters which indicates synthesizability of HEAs are also listed for each HEA. For the definition of δ, ΔH, Ω and VEC, see text and Refs. 1), 13), 33)–35).

Origin of the high Young’s modulus predicted for MoTcTaWRe in Table 1 can be traced back to its electronic structure. Roughly speaking the overall shape of the density of states (DOS) is determined by the lattice structure and the characteristics of BCC structure is large and deep valley structure appearing in the middle of the DOS. Basically, the valley structure distinguishes bonding states from anti-bonding states. Therefore, when the Fermi level is located at the valley the structure becomes more stable leading to large Young’s modulus. This is why transition elements nearly half-filled frequently appear in Table 1. Among them W, Re and Os are especially noticeable. In general, due to the larger delocalization of d-orbitals the hybridization becomes larger for later periods. This partly explain why W, Re and Os are frequently appear in Table 1. General trend between the calculated Young’s modulus and the DOS at the Fermi level in BCC HEAs was analyzed by Hayashi et al.¹³⁾

3.2 Atomic arrangement in HEAs at finite temperature

Table 2(a) summarizes the calculated EPIs in CrMnFeCoNi (so called Cantor alloy³⁶⁾). EPIs for all combinations of constituent elements are calculated up to the 3rd nearest neighbors. As shown in Table 2(a), the nearest neighbor (1st n.n.) EPIs involving Cr are larger compared to the others. Since the 1st n.n. EPIs with Cr are positive, Cr has a tendency to avoid 1st n.n. configuration and Cr is excluded from the 1st n.n. site of Cr. On the other hand, for the 2nd nearest neighbor (2nd n.n.) EPIs involving Cr are somewhat smaller compared to the others. This means Cr favors to locate at 2nd n.n. of Cr.

Table 2 Calculated EPIs for (a) CrMnFeCoNi and (b) CrMnFeCoCu in the unit of meV.

The MCS results for CrMnFeCoNi by using the EPIs in Table 2(a) are shown in Fig. 1(a)∼(c). As seen in Table 2(a), due to the large magnitude of the EPIs involving Cr, SRO for the pairs including Cr show peculiar temperature dependence, namely below 1000 K 1st n.n. SRO approaches to 1 (Fig. 1(a)) and 2nd n.n. SRO approaches to −4 (Fig. 1(b)). This behavior corresponds to the arrangement of Cr atoms at the corner of conventional BCC unit cell resulting in L1₂ like ordering. The similar ordering was reported in previous theoretical works.¹⁷^,¹⁹⁾ For the elements other than Cr, no significant order is observed except for Ni which shows clustering at very low temperature. These typical arrangement of Cr and Ni can be seen in the snapshot of the atomic arrangement in CrMnFeCoNi as shown Fig. 1(c). In the experiment, homogeneous CrMnFeCoNi was successfully fabricated after the annealing at 1273∼1423 K for 1∼3 days.¹⁶^,¹⁸⁾

Fig. 1

Temperature dependence of SRO of (a), (d) 1st n.n., and (b), (e) 2nd n.n. (c), (f) the snapshots of the atomic arrangements taken at 500 K for both systems. Upper panels show the results for CrMnFeCoNi and lower panels for CrMnFeCoCu.

The calculated EPIs of CrMnFeCoCu are summarized in Table 2(b). The EPIs involving Cr shows similar behavior to the ones in CrMnFeCoNi, namely large positive values for the 1st n.n. and negative small values for the 2nd n.n. with the exception of 1st n.n. Cr–Cu EPI. The sign of the EPI for Cr–Cu pair is opposite to the case in CrMnFeCoNi. It is seen that the EPIs involving Cu are all negative for the 1st n.n. This means that Cu has a tendency for the cluster formation. As a result, in CrMnFeCoCu, Cu clustering should be observed for lower temperature. Actually this tendency of Cu clustering is significant as shown in the simulation results of SRO (Figs. 1(d) and (e)) and snapshot of atomic arrangement (Fig. 1(f)). The simulated results well correspond to the experimental observations of Cu clustering in CrMnFeCoCu.¹⁶⁾

4. Conclusion

In this paper, we have demonstrated two CMD concerning to the HEAs. Firstly, in order to overcome combinatorial explosion of material space of HEAs we have constructed prediction model of elastic constants. B, c′, c₄₄ were calculated by using the FPKKR-CPA method for 2555 BCC HEAs, and by the linear regression, we constructed models for the 3 elastic constants. For accurate calculations of elastic constants, the use of FPKKR is mandatory. Moreover, the CPA considerably reduces and simplifies the calculations electronic structures and elastic properties of HEAs. The LIDG method arrows us to generate independent higher order descriptors starting from rather simple basic descriptors. By optimizing the choice of the descriptors, we can construct linear models with the accuracy comparable to the standard neural network method. By using the constructed models, we have proposed several new HEAs with high Young’s modulus.

Secondly, we have simulated the atomic arrangement in HEAs by performing the MCS. The HEAs are modeled by using the Potts-like model and the interaction parameters were determined by the GPM. The simulated atomic arrangements in CrMnFeCoNi and CrMnFeCoCu reasonably agree with the experimental observations.

Since the (FP)KKR-CPA method can effectively provide the accurate information of configuration averaged properties of disordered alloys, this method is suitable for HEA systems. Moreover, a combination with machine learning methods or Monte Carlo simulations, we can extend CMD for larger materials space and for finite temperature to contribute the acceleration of searching new functional materials.

Acknowledgments

The authors would like to acknowledge Dr. Masako Ogura for providing us with her FPKKR-CPA code, Dr. Hitoshi Fujii for his guidance to the usage of the LIDG package. This work is partly supported by JSPS KAKENHI (Grant No. 20K05303, 18H05212), JST CREST (Grant No. JP-MJCR18I2) and MEXT as “Program for Promoting Researches on the Supercomputer Fugaku” (DPMSD, Project ID: JPMXP1020200307), and “Data Creation and Utilization-Type Material Research and Development Project” (Digital Transformation Initiative Center for Magnetic Materials, Grant No. JPMXP1122715503).

REFERENCES

Corresponding author

Register with J-STAGE for free!