Prediction model of elastic constants of BCC high-entropy alloys based on first-principles calculations and machine learning techniques

ABSTRACT By combining first-principles electronic structure calculations and machine learning techniques, prediction models of elastic constants are constructed for BCC high-entropy alloys (HEA) containing 5 different elements chosen from 3d, 4d and 5d transition metals with equal concentration. Three independent elastic constants of randomly selected 2555 HEAs are calculated by using the full potential Korringa–Kohn–Rostoker (FPKKR) method with taking configurational disorder into account within the coherent potential approximation (CPA). From the obtained database of the elastic constants, prediction models are constructed by the linear regression using the descriptors generated by the linearly independent descriptor generation (LIDG) method. By optimizing the selection of descriptors based on the genetic algorithm (GA), prediction errors of 10.2 GPa, 4.5 GPa, 2.4 GPa and 7.7 GPa are achieved for bulk modulus , shear moduli , and Young’s modulus , respectively. By using the generated model we propose some HEAs with low . It is well known that the magnitude of is closely related to the shape of the calculated density of states (DOS). This statement is reconfirmed within the BCC HEAs, i.e., HEAs with larger DOS at the Fermi level shows smaller Young’s modulus and vice versa. Graphical abstract


Introduction
High-entropy alloys (HEAs) are a new class of metallic alloys composed of more than 5 different elements with equiatomic or near equiatomic composition. Following the original works of the fabrication of HEAs [1][2][3][4][5], significant efforts have been devoted to synthesizing HEAs and investigating their physical properties from various points of view. So far, in this new category of alloys many exceptional physical properties have been reported [6], such as high hardness and high-temperature strength [7], good temperature dependence of strength and ductility [8], thermal stability [9], high hardness and low-density refractory [10], resistances to wear, corrosion, oxidation, and fatigue [11][12][13][14], bio-compatibility [15], and even superconductivity [16]. Thus, HEAs become a new category of alloy studies and provide opportunities to explore new functional alloys.
The concept of the HEA is different from conventional alloys composed of one or two major elements with additional impurity elements to control their physical properties. The essence of HEAs is the disorder and combinatorial complexity. For example, the configurational entropy of n-component HEA is calculated as S ¼ k B log n, and S reaches 1.61k B for n ¼ 5 with a completely random configuration. Due to this large configurational entropy of the system, the formation of a random solid solution phase is expected, as long as the magnitude of mixing enthalpy is small enough. This allows us to choose any combinations of constituent elements for fabricating HEA and provides us with more freedom in the choice of the materials, for example, suppose 5 constituents are chosen randomly from 25 transition elements, there are 53,130 possibilities for HEAs. As for material properties, the randomness of the local atomic configuration is supposed to affect the dynamics of impurities and defects leading to favorable mechanical properties in HEAs. An expectation to realize nontrivial effects arising from a novel combination of elements is called the cocktail effect and this is why the HEA attracts much attention from materials scientists [3][4][5].
HEAs have numerous combinations of elements, so more HEAs with good properties will be discovered and useful for many applications in the future. However, randomness and combinatorial complexity inherent in HEAs make it difficult to establish a guideline for controlling the materials properties of HEAs. Since exhaustive experimental search is not practical, computational approaches have been employed to support experimental investigations on functional HEAs. Among the computational methods, first-principles electronic structure calculations have been one of the most reliable methods due to its ability to simulate materials properties by tracing back their electronic origin [17,18]. Most of the first-principles calculations were done based on the density functional theory (DFT) by using various band structure methods. So far, for the HEA researches, the DFT calculations were applied to thermodynamic properties [19][20][21][22][23][24][25], local distortions [26,27] phase stabilities [18,28], mechanical properties [29][30][31][32][33][34][35] and magnetism [36][37][38]. However, due to its computational costs, computational materials design based on the DFT calculations is limited for HEA systems and still at its development stage.
In this paper, as one of the fundamental mechanical properties, we focus on the elastic constants of equiatomic quinary HEAs in BCC structure and try to construct a prediction model with reasonable predictability so that we can make use of it for the computational materials design of HEAs. Our strategy is along with the idea of 'Materials Informatics' [39][40][41][42][43], namely, first we prepare a database of the elastic constants of HEAs based on the DFT calculations, then construct a prediction model of the elastic constants by applying machine learning techniques to the database [44]. For the machine learning, we use simple linear regression with descriptors generated by the linearly independent descriptor generation (LIDG) method developed by Fujii et al. [45,46]. By tuning the selection of descriptors with employing the genetic algorithm (GA), optimization of the model is carried out.
In the linear regression method, the objective variables are expressed by simple polynomials of the combinations of several descriptors. Therefore, the generated model can be interpretable. Taking this advantage of the linear model, we will discuss physical trends from the relation between the elastic constants and the selected descriptors. As for the machine learning technique, the neural network (NN) is widely used and well accepted as standard tool to realize machine learning. In spite of its usefulness for constructing accurate prediction model, it is sometimes criticized that the interpretation of the NN model is not straightforward. Thus the NN modeling and the LIDG-GA modeling are complementary with each other. In the present paper, the prediction model is also constructed by using the NN and its performance will be compared to the LIDG-GA model. We also try to demonstrate the applicability of our model for the prediction of HEAs with desired elastic constants, low Young's modulus (E) in the present work. Additionally, we will analyze the electronic structure of predicted low-E HEAs and the relation between electronic structure and low-E will be investigated.

Elastic constants
For the present procedure of the modeling, to construct a database for machine learning the elastic constants of HEAs should be calculated as many as possible. Due to the requirement for the computational cost for taking an average of physical properties over possible atomic configurations, it is impractical to apply the standard band-structure method for supercells of the special quasi-random structure (SQS). In this paper, for the calculations of the electronic structure and elastic constants of HEAs, we employ the fullpotential Korringa-Kohn-Rostoker (FPKKR) method with combining the coherent potential approximation (CPA). By using the FPKKR-CPA method, the configuration average of the electronic structure of HEAs can be calculated effectively from first-principles. For the present calculations, we use the program package developed by Ogura et al. [47]. Since the present implementation of the CPA to the FPKKR method depends on the single-site approximation, local fluctuation of atomic configurations and local lattice relaxation are not considered. It was pointed out that even for this simplification the CPA gives reasonable estimation for the fundamental properties of HEA including the lattice constant, relative structural stability between BCC and FCC and the elastic properties [17,28,32,37], and shown to be applicable for computational design of HEAs. On the other hand, it was also pointed out that the CPA gave a considerable error in the estimation of the heat of formation. The error reaches up to 30 meV/atom for MoNbTaVW, therefore the prediction of the formation of the solid solution by the CPA should be separately examined [18]. In the present paper, we do not focus on the modeling of the heat of formation, and in sec 3.2 we will use empirical rule for the screening of HEAs.
All of the present FPKKR-CPA calculations were spin-polarized and performed within the scalar relativistic approximation without the spin-orbit interaction. The exchange-correlation energy functional proposed by Perdew, Burke and Ernzerhof was employed [48]. The maximum angular momentum for the expansion of the Green's function is l max ¼ 4, and the 20 � 20 � 20 mesh points are prepared in the first Brillouin zone and the irreducible ones are used for the k-space integration. The electron density was calculated by taking the imaginary part of the Green's function integrated by using 35 energy mesh points on the complex energy contour. The width of the contour is set 1:0,1:5 Ry depending on the system. For a cubic structure, there are 3 independent elastic constants c 11 ; c 12 and c 44 . To determine these 3 constants, we apply 3 different kinds of deformations. The first one is the change of the unit cell volume V. By fitting the Murnaghan's equation of state (equation 1) to the calculated total energy E as a function of V, equilibrium volume V 0 , ground-state total energy E 0 , bulk modulus B and its derivative B 0 are determined. In the present calculations, for the determination of B volume change of À 3%, þ 3% around V 0 was assumed.
Then, we apply orthorhombic deformation δ o and monoclinic deformation δ m [31,32] with keeping V ¼ V 0 , namely, primitive translation vectors are transformed with the following matrix, respectively, The total energy change ΔE due to volumeconserving orthorhombic and monoclinic deformations is calculated as a function of δ o and δ m and fitted to the following equations 2 and 3 to derive two cubic shear moduli c 0 and c 44 , respectively. B and c 0 are related to the elastic constants as B ¼ ðc 11 þ 2c 12 Þ=3 and c 0 ¼ ðc 11 À c 12 Þ=2. In the present calculations, δ o and δ m are changed from 0,0:03 for the determination of c 0 and c 44 .
For constituent elements of HEAs, we assume 3d transition metals (Sc, Ti, V, Cr, Mn, Fe, Co, Ni and Cu), 4d transition metals (Y, Zr, Nb, Mo, Tc, Ru, Rh, Pd and Ag) and 5d transition metals (Ta, W, Re, Os, Ir, Pt and Au). By choosing 5 elements randomly from these 25 transition metals we generate in total 2555 HEAs and calculate electronic structure and derive V 0 ; B; c 0 and c 44 . Once the elastic constants are obtained Young's modulus E can be derived as where G is the share modulus. We assume poly-crystalline sample and G is supposed to be arithmetic mean of Voigt average G V and Reuss average G R , where G V ¼ 2c 0 þ3c 44 5 and G R ¼ 5c 0 c 44 2c 44 þ3c 0 . For typical HEAs, calculated Young's moduli are compared with the previous theoretical results and experimental values in Table 1, and it is confirmed that they show reasonable agreement with each other.

Construction of prediction model
For model construction, we employ general linear regression method. In this case, the choice of descriptors is of primary importance. Basically, we follow the procedure applied to sparse modeling of chemical bonding in binary compounds by Kanda et al. [46]. In the present study, objective variables are bulk modulus B, elastic constants c 0 and c 44 of a target HEA, and as for the basic descriptors we use the calculated physical values of constituent elements of the HEA, namely, B; c 0 ; c 44 , lattice constant a, group of element g, period of element p, atomic number Z and electron density parameter r s . r s is defined as 4π 3 r 3 s ¼ 1=ρ, where ρ is electron density in the interstitial region of the unit cell [53]. The electron density in the interstitial region is calculated from the number of electrons outside the muffin-tin (MT) spheres and the volume of the interstitial space between MT spheres. For the estimation of interstitial electron density, we performed KKR calculation within the MT approximation by using the Akai-KKR method [54,55] with using the lattice constant obtained by the present FPKKR calculations. The MT radius is set so that the MT spheres touch with each other. All of the basic descriptors except for r s summarized in Table 2 are calculated by the FPKKR with assuming BCC structure for all elements irrespective of their equilibrium structure. For some elements, BCC is not stable and this is why the calculated shear moduli become negative or very small for some cases.
From the basic descriptors, first-order descriptors which are a kind of seeds for constructing higherorder descriptors are generated. Since a model of a HEA should be symmetric under the permutation of constituent elements, we set up the following firstorder descriptors which have permutation symmetry. Namely, concerning the physical value X, its average hXi ¼ P i¼1;5 X i =5 and the standard deviation σ X ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi q and the inverse of the average 1 hXi , where X i is the physical value X of the i-th constituent element. In total there are 25 first-order descriptors with including a constant term.
Then we generate higher-order descriptors by constructing products among the first-order descriptors up to 2nd order. From the set of all possible products, linearly independent ones are picked up by using the LIDG method [45,46]. Finally, from 325 2nd order products, 316 descriptors are chosen as linearly independent ones and used for the following linear regression. The details of the LIDG algorithm and python package developed by Fujii et al. can be found in ref [56].
It is not trivial which ones out of 316 descriptors play an important role in the prediction of elastic constants. Moreover, more descriptors do not necessarily mean the better model. Therefore, in this work, a selection of descriptors based on the cross-validation is performed. However, it is impossible to perform an exhaustive search for the optimized model, since the number of combinations of N descriptors out of 316 candidates can be too large to handle. Instead, we start with randomly chosen N descriptors and update the combination of descriptors by using the GA to find the best one which gives the smallest validation error. As for the estimation of validation error, we employ the leave-one-out cross-validation and monitor the root mean square error (RMSE) R ¼ RðNÞ ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi is the predicted residual for the i-th HEA by the model constructed with excluding the data of the ith HEA. In the present implementation of GA, population size is 2000 at most and updated to minimize R. The best 1000 sets are conserved. From the rest, 999 sets are generated by crossover and 1 set by mutation. The probabilities of crossover and mutation are 90% and 100%, respectively. After 500 updates, we choose the one with the  smallest R as the best model with N descriptors. By performing this GA optimization for each N, finally, the optimized model is obtained. For the learning process, the database of the calculated elastic constants is divided into two groups, and 80% of the data are used as a training and validation set and 20% of the data are used as a test set. For the model construction explained above, only the training and validation set is used, and the test set is used only for the final assessment of the model after the model construction.

Optimization of prediction models
Carrying out the LIDG and GA procedure described in Sec.2, we perform model construction for elastic constants of HEAs. In Figure 1, the RMSE is plotted as a function of number of descriptors N, for bulk modulus B, elastic constants c 0 and c 44 . For each N the model is optimized by GA as explained in the previous section, and RMSE is evaluated by using leave-one-out cross validation for the training-validation data set. The model constructions are performed independently for these 3 physical values. It is found that for B, c 0 and c 44 , the RMSE becomes minimum at N ¼ 160; 208 and 198, respectively. In Table 3, the obtained smallest RMSE calculated by using test data that is not used in the model construction are summarized as a prediction error. For comparison, we also perform modeling by using a neural network (NN) [57]. For the NN modeling, the 25 first-order descriptors prepared for the LIDG procedure are used as inputs and predict B, c 0 and c 44 with setting one hidden layer which has 64 neurons and a constant term. As the activation function we employed rectified linear unit (ReLU) function. Obtained RMSE by the NN model are also summarized in the table. It is found that the present LIDG-GA method offers an accurate model whose prediction error is comparable to the one obtained by the NN.
To assess the performance of the present models visually, in Figures 2(a)~(c), B, c 0 and c 44 of HEAs included in the test data are plotted on the 2d-plane by using FPKKR-CPA values and LIDG-GA predictions. By assuming a poly-crystalline sample, Young's modulus E is calculated for each HEA whose B, c 0 and c 44 are predicted to be positive and also plotted in Figure 2(d). As shown in the figures, the performances of the present LIDG-GA models are reasonable and successfully predict elastic properties of HEAs. Interestingly, similar performance is observed even for derived physical value E (RMSE = 7.7 GPa).
Since the present method is basically the linear regression, the constructed model is merely a polynomial of descriptors with materials parameters as explained in Sec.2. Therefore, it might be possible to discover the physical trend by analyzing which terms contribute to the prediction significantly. This possibility of interpretation of the constructed model is one of the advantages of using LIDG-GA [46]. For this purpose, we focus on the  model which has a smaller number of descriptors. As recognized from Figure 1, for all three objective variables predictability is already reasonable for the models with 5 descriptors, and the RMSE are evaluated as 13.0, 13.5 and 4.7 GPa for B, c 0 and c 44 , respectively. In Table 4, we list up the 5 descriptors for the small models of B, c 0 and c 44 . Note that almost all of the selected 5 descriptors are 2nd order products. This means that if we take only the first order descriptors, the degree of freedom of the model is very limited and to achieve reasonable accuracy by describing complicated dependence of objective variables on the descriptors we need to increase the complexity of the model. The beneficial point of the LIDG-GA method is its balanced nature between accuracy and interpretability. We can keep the model as simple as possible, and at the same time, the complexity of the model can be systematically controlled by increasing the order of the products in the model.
For the constructed small models we discuss which descriptors are important for determining these physical values. In Table 4, the importance of each descriptor is also summarized. We define the importance of j-th descriptor x j in the regression of physical value y as  descriptor x j . Namely, G 2 j directly means how much the error (residuals) increases when regression is performed without x j . By definition, G 2 j � 1 and the descriptor which has larger G 2 j is more important [45]. In the table, the obtained regression coefficients b are also summarized. It is reasonable that important descriptors for b and c 0 are related to hBi and hc 0 i, namely the averages of b and c 0 of the constituent elements. However, hc 44 i is not on the list of descriptors for c 44 . Since hBi has a large correlation between hc 44 i, the effect of hc 44 i might be included via hBi. This speculation is partly justified by the appearance of hBi in the list of descriptors for c 44 .
Since most of the selected important descriptors are the 2nd order products, direct interpretation of the model seems too complicated. Instead, we try to analyze Pearson's cross-correlation values between objective variables and the first-order descriptors which appear as the factors in the important descriptors listed in Table 4. As shown in Table 5, it is found that B correlates closely not only with hBi but also with the descriptors related to the lattice constants. Particularly, σ a shows a significant correlation with B and it has a negative regression coefficient. This means that the difference in the lattice constants of constituent elements becomes larger, the bulk modulus of the HEA becomes smaller. Such a negative correlation was also pointed out by Koval et al. [35].
In the case of c 44 , among the first-order descriptors, it is interesting to find a large correlation with hgi, the average of the number of groups of the constituent elements. In the present calculations the crystal structure is fixed to BCC, therefore, if we assume rigid band like behavior in the calculated DOS, the electronic structure is determined mainly by the average valence electron number which corresponds to hgi. This might be a reason why this descriptor controls the mechanical stability of HEA systems. However, in the case of c 0 , it shows a correlation with the inverse of hgi, thus it seems not to be straightforward to extract the physical origin of the relation between selected descriptors and elastic constants. By directly inspecting the DOS of HEAs the relation between mechanical stability and the electronic structure may become more clear. General discussion from this viewpoint will be presented in the next subsection.

Prediction of HEAs with low Young's modulus
Now, we try to demonstrate the application of the most predictable model constructed by the LIDG-GA to the materials design of HEA. As a target property, we focus on Young's modulus E of HEA and try to propose HEAs with a small E value. Note that the present models concern only elastic constants and do not predict the stability of the BCC random solution phase. To propose realistic candidates, we perform further screening by using empirical rule for the fabrication of HEAs [3,20,21], namely, δ < 8:5ð%Þ; À 22 < ΔH < 7ðkJ=molÞ; Ω > 1 and VEC < 6:87, where δ ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi P i¼1 c i ð1 À r i hri Þ 2 q � 100,  Table 6 shows the candidates of low-E HEAs predicted by the present model and screened by the above criteria for stability. Among them, TiVZrNbTa and TiZrNbMoTa [15] were already studied by experiments to fabricate bio-compatible materials, but their E were not reported yet.
To investigate the origin of low Young's modulus, the calculated electronic structure is analyzed for typical cases. For this purpose, we pick up two typical HEAs from the presently generated database of HEAs by the FPKKR-CPA calculations. One is VCrMoWRe which has the highest E of E ¼ 321 (GPa) and the other one is ScTiVZrNb which has the lowest E of E ¼ 46 (GPa). B, c 0 and c 44 of these two HEAs are all positive, and the predicted E of VCrMoWRe and ScTiVZrNb by the model are 308 and 45 (GPa), respectively.
As shown in Figure 3(a), for VCrMoWRe (high-E case), the DOS shows a clear dip and that the Fermi level E F is located at the dip. This characteristic dip structure in the DOS is typical for BCC structure and the dip distinguishes the bonding states from the antibonding ones. Thus, it is intuitively recognized that the high-E of VCrMoWRe might be related to the electronic stability of this HEA. This speculation is partly confirmed by looking at the DOS of low-E HEA ScTiVZrNb where the DOS shows peak structure at the E F as shown in Figure 3(b).
To investigate the relation between Young's modulus and the electronic structure, the correlation of the DOS at E F with E is plotted in Figure 4(a) for the 1590 HEAs predicted to have B > 0; c 0 > 0 and c 44 > 0 by the FPKKR-CPA. As recognized from the figure, the correlation between the DOS at E F and E is easy to recognize, and the calculated Pearson correlation coefficient is -0.48. A similar negative correlation with E F is also found with c 0 for the 2555 HEAs calculated by the FPKKR-CPA in the database (Figure 4(b), the correlation coefficient of -0.51). The correlation between the DOS at E F and structural instability is well known and the martensitic transformation of Fe-Pt alloys were investigated from this view point [58]. Therefore, the present finding concerning to the Young's modulus is not entirely new idea. It can be regarded as a new justification of that idea within the new class of metallic alloys BCC HEA with a wide range of combinations of constituent elements.
This information concerning to the correlation between the DOS at E F and E might be useful to construct a new prediction model for the mechanical properties of HEAs. Of course, to use such a model, in which the electronic information (DOS at E F ) is included in the set of descriptors, we need to perform a first-principles calculation for a target HEA. This might require a certain effort before obtaining

Summary
In this paper, we have demonstrated the construction of a prediction model by using the LIDG-GA method based on the first-principles database of elastic constants of HEAs in BCC structure. Our database was originally generated by performing the FPKKR-CPA calculations for 2555 HEAs whose constituent elements are randomly chosen from 3d, 4d and 5d elements. We have shown that even with the simple linear model we can achieve accurate models for B, c 0 and c 44 by choosing descriptors carefully. In the present study, the descriptors were generated from basic descriptors including elastic properties of B; c 0 ; c 44 , lattice constant a, group of element g, period of element p, atomic number Z and electron density parameter r s calculated by the FPKKR for each constituent element of the HEA. From these basic descriptors, owing to the LIDG method, higher-order descriptors are automatically generated with keeping their linear independence. We can optimize the combination of descriptors by using the GA to construct the most predictable model. In the present study, prediction errors measured by the RMSE are 10.2, 4.5 and 2.4 GPa for B, c 0 and c 44 , respectively. By analyzing the small model with 5 descriptors, we recognized the physical trend of B. Namely, HEAs with the smaller difference in the lattice constants of constituent elements show the larger B. The most predictable model was applied for the prediction of low Young's modulus HEAs. Experimental verification might be desirable for the present predictions. By investigating the electronic structure of proposed HEAs, the correlation between Young's modulus of HEA and the density of states at the Fermi level is figured out. This finding might be useful to construct another model of the physical properties of HEAs. The presently proposed computational method, namely the combination of FPKKR-CPA calculations and LIDG-GA modeling, might be one of the effective tools towards the computational design of functional HEAs.