RESEARCH ARTICLE | JUNE 07 2022 Toward machine learning for microscopic mechanisms: A formula search for crystal structure stability based on atomic properties  Udaykumar Gajera  ; Loriano Storchi ; Danila Amoroso ; Francesco Delodovici; Silvia Picozzi J. Appl. Phys. 131, 215703 (2022) https://doi.org/10.1063/5.0088177 CrossMark   View Export Online Citation 11 December 2023 10:28:07 Journal of ARTICLE scitation.org/journal/jap Applied Physics Toward machine learning for microscopic mechanisms: A formula search for crystal structure stability based on atomic properties Cite as: J. Appl. Phys. 131, 215703 (2022); doi: 10.1063/5.0088177 Submitted: 14 February 2022 · Accepted: 16 May 2022 · View Online Export Citation CrossMark Published Online: 7 June 2022 Udaykumar Gajera,1,2,a) Loriano Storchi,3,b) Danila Amoroso,1,4,c) Francesco Delodovici,1,5 and Silvia Picozzi1 AFFILIATIONS 1 Consiglio Nazionale delle Ricerche, CNR-SPIN c/o Università “G. D’Annunzio,” 66100 Chieti, Italy 2 Chemistry Department, University of Turin, via Pietro Giuria, 7, 10125 Torino, Italy 3 Dipartimento di Farmacia, Università degli Studi G. D’Annunzio, 66100 Chieti, Italy 4 NanoMat/Q-mat/CESAM, Université de Liège, B-4000 Liège, Belgium 5 Université Paris-Saclay, CentraleSupélec, CNRS, Laboratoire SPMS, 91190 Gif-sur-Yvette, France a) Author to whom correspondence should be addressed: uday.gajera@edu.unito.it b) loriano@storchi.org. URL: https://www.storchi.org/ c) danila.amoroso@uliege.be 11 December 2023 10:28:07 ABSTRACT Machine-learning techniques are revolutionizing the way to perform efficient materials modeling. We here propose a combinatorial machine-learning approach to obtain physical formulas based on simple and easily accessible ingredients, such as atomic properties. The latter are used to build materials features that are finally employed, through linear regression, to predict the energetic stability of semicon- ducting binary compounds with respect to zinc blende and rocksalt crystal structures. The adopted models are trained using a dataset built from first-principles calculations. Our results show that already one-dimensional (1D) formulas well describe the energetics; a simple grid- search optimization of the automatically obtained 1D-formulas enhances the prediction performance at a very small computational cost. In addition, our approach allows one to highlight the role of the different atomic properties involved in the formulas. The computed formulas clearly indicate that “spatial” atomic properties (i.e., radii indicating maximum probability densities for s, p, d electronic shells) drive the sta- bilization of one crystal structure with respect to the other, suggesting the major relevance of the radius associated with the p-shell of the cation species. Published under an exclusive license by AIP Publishing. https://doi.org/10.1063/5.0088177 I. INTRODUCTION Aflowlib,15 C2DB,16,17 QPOD,18 Materials Project,19 Materials Cloud,20 and related AiiDA,21,22 provide researchers with a huge collection of Modeling material properties with high accuracy and low compu- basic first-principles results. A large amount of ab initio data is thus tational cost is one of the grand-challenges in materials science and available, which can be used for deeper analyses and studies, provided engineering. The development of ab initio methods has provided accu- one can count on proper tools to extract relevant information out of rate tools for material properties prediction and their further optimiza- them. tion; nevertheless, one disadvantage of approaches relying only on In the last few years, materials scientists have developed differ- first-principles simulations is the high cost required in terms of com- ent machine-learning (ML) methods to rationalize the data putational resources and simulation time. In recent years, the continu- analysis.23–33 Each method has its own specific advantages and lim- ous growth of available computational power1 has stimulated scientists itations. Methods, such as random forest34 or neural network to move in the direction of high-throughput simulations.2–10 Along (NN),35 which is mainly behind the Deep Learning (DL), are very this line, open access databases, such as OQMD,11,12 NOMAD,13,14 efficient36 but not always transparent, partially blurring the J. Appl. Phys. 131, 215703 (2022); doi: 10.1063/5.0088177 131, 215703-1 Published under an exclusive license by AIP Publishing Journal of ARTICLE scitation.org/journal/jap Applied Physics comprehension of the role played by the input variables in the final discussed in what follows. Let us also remark that large packages results. Nonetheless, over the last few decades, improvements are nowadays available to the scientific community, already provid- toward the interpretability of such “black-box” ML models have been ing advanced and well-tested features for use in ML applied to made through additional methodologies,37 such as model-agnostic materials science.57,58 However, we here strictly followed the spirit methods, which in turn are divided into global and local interpreta- of Ghiringhelli et al.55 and, therefore, chose their same atomic fea- tion techniques (see Ref. 38 and references therein). For instance, out tures, as detailed below. of various global methods, we can cite: the permutation feature importance,39,40 which associates each feature importance values II. METHODOLOGY depending on how much the model error increases when its values are shuffled; the functional decomposition,41 which decomposes the The approach we present here can be regarded as a combina- complex prediction function into smaller parts; and the global surro- torial machine-learning: a set of basic atomic properties (APs, gate,42 which replace the original model with a simpler model that listed in Table S2 of the supplementary material) are randomly can be more easily interpreted. On the other hand, among the local combined (though under certain initial constraints detailed below) methods, we can cite: the local surrogate models (LIME),43 which to build a set of material features (MFs). The generated features are replaces the complex model with a locally interpretable surrogate then used to train a LR model, where the energy difference model, and the SHapley Additive exPlanations (SHAP),44 which is between rocksalt and zinc blende structures is the dependent vari- based on Shapley values and computes the contribution of each able (i.e., the label). Then, we select the best-performing model feature to the prediction. Also, in the specific case of DL, which struc- according to standard performance metrics, such as the root mean tures algorithms in multiple layers to create “artificial neural net- squared error (RMSE). The final result of this procedure is a works,” thus enhancing the complexity in the prediction’s “formula,” which is a concise and clear representation of the rela- interpretation, other specific interpretation methods have been tionship between the used atomic properties and the energy differ- proposed,45–47 in addition to the already cited model-agnostic ones. ence between RS and ZB phases. In the following, we describe in However, when targeted case studies allow, the easiest way to detail the different steps of our approach. achieve a deeper understanding of machine-learning results is to rely on interpretable models, such as linear regression (LR),48–51 logistic A. Dataset preparation and materials regression,52 and decision trees.53 This can apply here to our target As mentioned, we aim at predicting the total energy difference case, i.e., the prediction of the difference in total energy (ΔE) between (ΔE ¼ ERS  EZB ) between RS and ZB phases of cubic crystal struc- 11 December 2023 10:28:07 rocksalt (RS) and zinc blende (ZB) crystal structures in semiconduc- tures for 82 semiconductor binary AB compounds (the dataset is tor binary AB compounds, a prototypical case in materials science. reported in Table S2 of the supplementary material). We employed Accordingly, being our goal the creation of formulas linking the total energies reported in Ref. 55, which were calculated through target label (ΔE) to a set of basic atomic features, we selected the LR density functional theory (DFT)59,60 within the local density algorithm both in its one-dimensional and multi-dimensional forms. approximation (LDA61). In closer detail, we here propose a ML-based approach to build The construction of the material features is based on primary sets of features (or descriptors) starting from a given set of basic vari- atomic properties of the constituents, also taken from Ref. 55. To ables (e.g., atomic properties), which are subsequently used to con- facilitate the physical interpretation of each MF, the APs are sub- struct LR models (or formulas). The final outcome of our procedure is divided into two different kinds: (i) “energy” properties, including a transparent formula, not necessarily of easy mathematical formula- the highest occupied Kohn–Sham level (HOMO), the lowest tion, but revealing which part of the input mostly affects the output,54 unoccupied Kohn–Sham level (LUMO), ionization potential (IP), i.e., allowing the identification of the main driving physical features. electron affinity (EA) and (ii) “spatial” properties, including rs , rp , To test our method, we target a prototypical case in materials and rd , i.e., the radii where the radial probability density of the science: indeed, inspired by the original work of Ghiringhelli et al.,55 valence s, p, and d orbitals, respectively, reaches its maximum.4 we optimized our models to predict the difference in energy between RS and ZB; from that optimization, a classification of the most stable crystal structure between RS and ZB for semiconductor AB binary B. Formula construction compounds naturally derives. We rely on the LR48,49 approach to obtain a direct interpretation To identify useful features, we generate combinations of basic of the dependent and independent variables. The construction of a atomic properties (i.e., the independent variables in our approach) useful LR model can become troublesome, requiring a linear depen- of the material constituents through a combinatorial approach.56 dence between features. In Ref. 55, the authors implemented an auto- We then carry out an analysis of the emerging best-performing for- mated feature selection method employing the LASSO regression mulas, identifying the role of specific atomic features in determin- analysis method.55,63 In our work, we use a combinatorial approach to ing the final stabilization of the crystal structure. Finally, we test the generate the dependent variable (material features) to be used within predictive capability of the obtained formulas by applying them to the linear equations and thus to finally obtain the formulas. “new” compounds (i.e., outside the dataset used for training the In Fig. 1, we illustrate the workflow of the formula generation model), finding an overall satisfactory agreement with first- and selection using LR. The process starts with the selection of the principles results. As already mentioned, our approach is similar to APs to be combined. Afterward, we choose prototype functions what was originally proposed by Ghiringhelli et al.,55 though with that are simple analytical operations applied to the APs. In our some differences and further extensions, which will be carefully case, we selected five prototype functions, f (x), namely, J. Appl. Phys. 131, 215703 (2022); doi: 10.1063/5.0088177 131, 215703-2 Published under an exclusive license by AIP Publishing Journal of ARTICLE scitation.org/journal/jap Applied Physics FIG. 1. (a) Basic atomic properties (APs) used to construct the material features. (b) Crystal structures of RS and ZB (plot made using the VESTA tool).62 Gray (yellow) spheres represent A (B) atoms. (c) Workflow for formula construction, machine-learning methodology, validation, and MF selection procedures. In the AB compounds, A is the atom with the lowest electronegativity. pffiffiffi 11 December 2023 10:28:07 x, x2 , x3 , x, ex , where x is an AP. Then, we obtain the final set of • GEN4: combine two prototype functions with the same physical MFs by combining different prototype functions via the combina- dimensions at both the numerator and denominator, torial approach (see, for instance, Ref. 56) and applying the follow- ing additional set of rules: f1 (AP1 ) w f2 (AP2 ) MF ¼ , (4) f3 (AP3 ) w f4 (AP4 ) • GEN1: combine two prototype functions in the numerator, forcing them to belong to the same kind of APs, which is both where w¼ þ   4. “spatial”-like or both “energy”-like; one prototype function is at the denominator with the only constraint to be non-zero, such as Each one of these set of rules corresponds to a different MF generator. f1 (AP1 ) + f2 (AP2 ) From the implementation point of view, each generator is a MF ¼ : (1) Python64 function that produces a set of strings. Therefore, we can f3 (AP3 ) easily exploit the Python capability to parse a source code and run a Python expression (code) within a program65 to compute all the MFs’ • GEN2: combine two prototype functions with the same kind of values starting from the generated sets of strings. This allows for an APs at the numerator and a single prototype function at the easy implementation and plugin of other generators, as well as to denominator with an argument of a different kind with respect easily adopt different sets of atomic properties, leaving the workflow to the numerator ones. For instance, if AP1 in f1 (AP1 ) and AP2 unchanged: a new generator can be introduced implementing a in f2 (AP2 ) are “energy” terms (i.e., EA or HOMO), then AP3 Python function returning a list of strings, each one being a valid MF. must be a “spatial” term (i.e., rp ), Finally, in order to choose the optimal formula, we build a LR model for each of the generated MF. To practically select the best f1 (AP1 ) + f2 (AP2 ) model, i.e., the “best formula,” we randomly split the full dataset MF ¼ : (2) f3 (AP3 ) into 90% as a training set to train/initialize the model and 10% as a test set to check model’s performance. We perform this random • GEN3: combine two prototype functions at both the numerator splitting N times (with N ¼ 150) for each model, and we calculate and denominator without any constraints, the RMSE from the test set for each run. Afterward, we again verify the top ten resulting best formulas with a higher value of training set and test set splitting, with N ¼ 1000. We average it over all N f1 (AP1 ) + f2 (AP2 ) MF ¼ : (3) splitting, and we obtain avg(RMSE), as reported in the tables of f3 (AP3 ) + f4 (AP4 ) Sec. III. J. Appl. Phys. 131, 215703 (2022); doi: 10.1063/5.0088177 131, 215703-3 Published under an exclusive license by AIP Publishing Journal of ARTICLE scitation.org/journal/jap Applied Physics We mention that different metrics for evaluating regression f1 (AP1 ), f2 (AP2 ), f3 (AP3 ), f4 (AP4 ) are the prototype functions build models can lead to different formula ranking. In this work, we rank on the primary atomic properties APi ; and m and q are the slope the obtained models based on the lowest avg(RMSE) for direct or angular coefficient and intercept, respectively, recursively comparison with a previous work.55 determined upon LR. In Table II, we report the optimized, best-performing formula from the different generators; the top 10 formulas are reported in C. Formula optimization Table S1 of the supplementary material. In order to further improve the performance of our models, To benchmark our grid search, we also used automated we introduce an additional step, which we refer to as “formula coefficient-optimizing methods: Nelder–Mead,68 conjugate gradi- optimization.” Specifically, we focus on the top ten formulas ent (CG),69 Broyden–Fletcher–Goldfarbnno (BFGS),70 and trun- obtained using each generator and the subsequent LR, as described cated Newton method (TNC).71 Although the resulting sets of in Sec. II B. After that, we use a grid search to find the relative coefficients are different in terms of single values with respect to weights of each prototype function of the atomic properties [i.e., those obtained via the grid search, the ratios between them are each fi (APi )] within the formula. A first grid search ranging almost preserved as well as the associated RMSE. In particular, for between 1 and 1 with the increasing step of 0.1 is used simultane- the case of GEN1 and GEN2, the ratio between the numerator coeffi- ously for all the weight coefficients (i.e., an exhaustive search cients a and b is preserved, and for GEN3 and GEN4, also, the through the specified subset of values for a, b, c, d coefficients is denominator coefficient ratio, between c and d, is preserved. In simultaneously performed). We multiply each fi (APi ) of the formula Fig. S3 of the supplementary material, we show the evolution of the by the weight coefficient, and we optimize the final RMSE value. Once RMSE and different ratios for different methods using 1D feature the procedure finds a set of optimal weight coefficients, two subsequent generated by GEN3. grid searches, with reduced incremental step values (0.01 and 0.001, Finally, we would like to underline that the whole proce- respectively) and a range of search are performed to obtain the final set dure, i.e., formula construction and its optimization, is not too of refined weight coefficients. Of note, for each set of weight coefficients expensive from a computational point of view. Indeed, as generated during the grid search, we also run the linear regression. reported in Table S4 of the supplementary material, for almost Thus, we are performing a proper formula optimization, as at each step all the generators, the whole computation can be performed in of the grid search, we are updating both the weight coefficients as well less than 4 h on a standard PC. Only GEN3, where 1 091 200 as the slope and intercept coming from the LR. In addition, it is impor- different formulas are generated and evaluated, is more time 11 December 2023 10:28:07 tant to underline that we made sure to find the global minimum solu- consuming (i.e., almost 15 h are needed). However, we under- tion when analyzing the N–D maps over the phase space of the line that the 1D formula construction procedure can be easily parameters (in the specific case of GEN2, we also double-checked the parallelized in order to drastically reduce its computational optimization results through an analytical minimization approach via burden. Mathematica66,67). To further clarify the procedure, we show here an exemplary equation, D. Higher-dimensional features a  f1 (AP1 ) w b  f2 (AP2 ) Following the idea of Ghiringhelli et al. (see Ref. 55), we also ΔE ¼ m  þ q, (5) c  f3 (AP3 ) w d  f4 (AP4 ) build 2D and 3D formulas as follows: we combined in all possible ways two or three different 1D MFs extracted from the best 1000 where ΔE is the targeted material feature (MF); a, b, c, d denote the ones and checked the avg(RMSE) using multiple LR for N test- weight coefficients scanned during the grid search; train set splits. Thus, the final equations, which relate ΔE to the TABLE I. 1D formulas, along with related statistics: avg(RMSE) denotes the root mean squared error for average over 1000 random train-test splits of dataset. Instead, the RMSE is the root mean squared error for the entire dataset as training and test. Similarly, the R2 values are calculated considering the entire dataset, and they show the quality of fit between predicted and actual values. The success rate (in percent) shows how many RS or ZB phases out of 82 have been correctly identified by the descriptor. The “Generator type” column indicates the different generators used to produce the corresponding descriptor. RMSEs are in eV. Formula avg (RMSE) RMSE R2 Success rate (%) Generator type 0:117  EA(B)IP(B) rp (A)2  0:342 0.1455 0.1423 0.89 89 1D descriptor55 r (B)3 exp[r (B)] 0:751  p r (A)2 s  0:317 0.1296 0.1193 0.92 90 GEN1 pffiffiffiffiffiffiffiffiffiffip pffiffiffiffiffiffiffiffiffiffiffi jIP(B)jþ jEA(A)j 0:285  rp (A)2  0:387 0.1367 0.1309 0.91 91 GEN2 p ffiffiffiffiffiffiffiffiffiffi rp (B)þ jrd (A)j 0:774  r (A)3 þr (B)3  0:303 0.0995 0.0963 0.95 94 GEN3 p p 1:155  rs (B)þrs (A) rp (B)3 þrp (A)3  0:368 0.1103 0.1058 0.94 96 GEN4 J. Appl. Phys. 131, 215703 (2022); doi: 10.1063/5.0088177 131, 215703-4 Published under an exclusive license by AIP Publishing Journal of ARTICLE scitation.org/journal/jap Applied Physics TABLE II. 1D formulas after the optimization step, along with related statistics. Notation as in Table I. RMSEs are in eV. Formula avg (RMSE) RMSE R2 Success rate (%) Generator type 0:127  0:800EA(B)1:000IP(B) 1:110rp (A)2  0:352 0.1457 0.1419 0.89 89 1D descriptor55 pffiffiffiffiffiffiffi 0:801 rp (B)0:606exp[rp (A)] 1:870  1:010rp (A) 3  0:968 0.1191 0.1143 0.93 91 GEN1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0:876 jHOMO(B)jþ0:468 jLUMO(B)j 0:477  1:110rp (A)2  0:372 0.1340 0.1296 0.91 91 GEN2 pffiffiffiffiffiffiffiffiffiffi 0:642rp (B)þ0:502 jrd (A)j 1:609  1:170r (A)3 þ1:170r (B)3  0:309 0.0991 0.0961 0.95 94 GEN3 p p 0:878rs (B)þ0:200rp (A) 1:207  0:512rp (B)3 þ0:610rp (A)3  0:359 0.1045 0.1016 0.94 99 GEN4 basic atomic features, are written as follows: comparing values obtained from the ML-predicted ΔE formula with the corresponding ab initio calculated values. In closer detail, a1  f1 (AP1 ) w b1  f2 (AP2 ) we focused on different alloys, obtained by changing, respectively, ΔE ¼ m1  c1  f3 (AP3 ) w d1  f4 (AP4 ) the concentration of A-site atoms, such as [Ax A01x ]B, and of B-site a2  f5 (AP5 ) w b2  f6 (AP6 ) atoms, such as A[Bx B01x ]. Accordingly, one can test the efficiency þ m2  þq (6) of the formulas by checking the energy difference for intermediate c2  f7 (AP7 ) w d2  f8 (AP8 ) concentrations as obtained from optimized 1D formulas and for the general 2D formulas and compare their trend with respect to first-principles results. To this end, ab initio electronic-structure simulations were carried out a1  f1 (AP1 ) w b1  f2 (AP2 ) within DFT and LDA functional. Calculations were performed ΔE ¼ m1  using the VASP72–74 code, employing a 8  8  8 k-mesh for the c1  f3 (AP3 ) w d1  f4 (AP4 ) Brillouin zone sampling. We verified that the results obtained with a2  f5 (AP5 ) w b2  f6 (AP6 ) þ m2  the pseudopotential VASP for the parent binary compounds were 11 December 2023 10:28:07 c2  f7 (AP7 ) w d2  f8 (AP8 ) consistent with those reported by Ghiringhelli et al., calculated a3  f9 (AP9 ) w b3  f10 (AP10 ) with the all-electron FHI-aims code.75 For simulations at different þ m3  þq (7) c3  f11 (AP11 ) w d3  f12 (AP12 ) concentrations, we adopted the so-called virtual crystal approxima- tion (VCA) based on virtual atoms interpolating between the real for the 3D ones. The comparison between performance is discussed constituent atoms.76,77 However, as well known from the literature, in Sec. III. the VCA approach neglects some effects, such as local distortions around atoms and, as such, should not be expected to reproduce fine details of disordered alloy properties.78 Accordingly, in some E. Test of predictive power of the ΔE formula for novel cases (i.e., for Mgx Ca1x Se alloys), in order to mimic disordered AB compounds structures with an improved accuracy, we calculated total energies After obtaining the optimized 1D formulas for ΔE in the case using supercell structures, rather than using the VCA method on of AB compounds, we aimed at further verifying their validity and primitive unit cells. Specifically, the considered supercell is the predictive power by considering additional AB systems (i.e., which cubic unit cell composed of four AB formula units with planes of were not originally included in the ML training set) and by cations alternating along the c direction (see Fig. S4 in the TABLE III. 2D formulas, along with related statistics. Notation as in Table I. RMSEs are in eV. Formula avg (RMSE) RMSE R2 Success rate (%) Generator type jrs (A)rp (B)j 0:113  EA(B)IP(B) rp (A)2  1:558  exp[rs (A)]  0:133 0.1041 0.0988 0.95 96 2D descriptor55 pffiffiffiffiffiffiffiffiffiffi rp (B)3 exp[rp (A)] rp (A)  jrd (A)j 2 0:342  rp (A)3  1:042  exp[rp (A)]  0:062 0.0989 0.0944 0.95 89 GEN1 pffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi IP(B)þ jIP(A)j r (A)  3 r (A) 0:081  rp (A)3  0:001  exp(HOMOKS(A)) s d  0:062 0.1163 0.1100 0.93 86 GEN2 pffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffi r (A) jrd (A)j rs (B)þ jrp (B)j 1:175  pr (B)3 þr (A) 3 þ 0:513  rp (B)3 þrs (A)3  0:250 0.0911 0.0878 0.96 87 GEN3 s p pffiffiffiffiffiffiffiffiffi r (A) jr (B)j 0:618  d 3 pp ffiffiffiffiffiffiffiffi þ 1:097  r (B)3 þr (A)3  0:384 r (A)=r (B) p p 0.0995 0.0955 0.95 92 GEN4 rp (A)  rd (A) p p J. Appl. Phys. 131, 215703 (2022); doi: 10.1063/5.0088177 131, 215703-5 Published under an exclusive license by AIP Publishing Journal of ARTICLE scitation.org/journal/jap Applied Physics 3D descriptor55 Generator type GEN1 GEN2 GEN3 GEN4 Success rate (%) 93 90 91 98 93 0.97 0.95 0.92 0.96 0.96 R2 0.0756 0.0933 0.1205 0.0834 0.0919 RMSE avg (RMSE) 0.0818 0.1003 0.1300 0.0875 0.0989 FIG. 2. Energy difference between rock salt and zinc blende, ΔE (in eV), as a function of rp (A) for different binary compounds (blue triangles). Data fit func- tions are also shown, using proportionality to rp (A)2 and rp (A)3 via a green dashed line and a red straight line, respectively. The inset shows the linear fit of ΔE for f (rp (A)) [ [0 : 1] as a function of 1=rp (A)2 (green triangles) and 1=rp (A)3  0:310 (red dots). Resulting slope (m) and intercept (q) are, respectively, 1.120 (eV Å2 ), 0.360 (eV) for 1=rp (A)2 and 1.184 (eV Å3 ), 0.240 (eV) for 1=rp (A)3.  0:016  HOMO(A)exp(LUMO(B)) 11 December 2023 10:28:07 (rp (A)3 ) exp(rp (A))exp(rd (B))  0:358 0:885  r (A)2 þr (A)3  0:417  r (A)3 þr (B)3  0:579  r (B)2 þr (A)3  0:616  0:124  p r (A)3 d  1:87 TABLE IV. 3D formulas, along with related statistics. Notation as in Table I. RMSEs are in eV. pffiffiffiffiffiffiffiffiffiffi rp (A) jrd (A)j pffiffiffiffiffiffiffiffiffiffi r (B)2  jr (A)j IP(A)2 EA(A)2  0:023 s p p jrp (B)rs (B)j exp[rd (A)] jEA(B)jþexp(EA(B)) þ0:038   3:782  rs (A)3 pffiffiffiffiffiffiffiffiffiffiffi rs (A)exp[rs (B)] pffiffiffiffiffiffiffiffiffi r (A)  jr (B)j ffi r (B) jrd (A)j pffiffiffiffiffiffiffiffiffi 0:635  r (A)3 þr (B)3 þ0:730  pr (A)3 þr (B) 3 p 0556  p r (A)3 p þ 0:364  p exp[rp (A)]d p jrs (A)rp (B)j þ 0:266  exp[rs (A)] 2 s p  1:806  (LUMO(A)þHOMO(B)) rp (B)exp[rp (A)] ffi pffiffiffiffiffiffiffiffi IP(B)= IP(A) rp (A)3 r (B) exp[r (A)] FIG. 3. Total energy difference ΔE as a function of concentration (x) for p EA(B)IP(B) p [Cax Sr1x ]Se, [Mgx Ca1x ]Se, and [Bex Mg1x ]Se alloys, highlighted in blue, rp (A)2 pffiffiffiffiffiffiffi p green, and pink regions, respectively. Energy differences are predicted using p 3 original and optimized 1D descriptors constructed using GEN3 and GEN4 and p 0:056  Formula verified using DFT (black line with diamond points) within VCA. For an improved 0:108  accuracy, the two asterisk-highlighted intermediate points in the [Mgx Ca1x ]Se region are calculated using the supercell approach rather than VCA. J. Appl. Phys. 131, 215703 (2022); doi: 10.1063/5.0088177 131, 215703-6 Published under an exclusive license by AIP Publishing Journal of ARTICLE scitation.org/journal/jap Applied Physics FIG. 4. Total energy difference ΔE as a function of concentration (x) for Sr[Sx Se1x ], see panel (a), and B[Px Sb1x ] alloys, see panel (b), predicted from original and opti- mized 1D descriptors constructed using GEN3 and GEN4. Model predictions are verified using energy differences calculated via DFT59,60 (black line with diamond points). supplementary material). The k-mesh was modified accordingly to the rocksalt or the zinc blende structure, we extracted the best ten maintain the same density of points employed in the simulations of formulas with the lowest avg(RMSE) from each generator (so-called primitive cells. original formulas) and then apply the formula optimization, as 11 December 2023 10:28:07 detailed in Sec. II C. This procedure attributes relative weights to each f (AP), allowing us to measure the importance of the individ- III. RESULTS AND DISCUSSION ual atomic properties in driving the energy stabilization. In princi- In this section, we will analyze the final formulas as obtained ple, the avg(RMSE) value depends on random test-train splits that from different generators. The results are shown in Tables I–IV; in we perform to our dataset. Therefore, to reduce the effect of ran- the first row, we report the results obtained by Ghiringhelli et al.55 domization, as target model performance metrics, we rank our for comparison. optimized formula based on the RMSE of the whole dataset, rather First, by comparing the avg(RMSE) values, we note that all 1D than based on avg(RMSE). By comparing Tables I and II, it is formulas obtained from our different generators better perform evident that the optimization procedure can further change the with respect to the 1D ones reported in Ref. 55, where the authors formula ranking, providing a different final “best formula” with used the automated feature selection method LASSO.63 respect to the non-optimized formulas. In particular, we notice an Remarkably, some atomic primary features appearing in 1D formu- improvement in RMSE around 5%–10% after the formula las of Ref. 55 also appear in our obtained list of 1D formulas using optimization. GEN1 and GEN2; nevertheless, those are characterized by a higher Interestingly, our results reveal the size of the A cation to play avg(RMSE) than other formulas we obtained via our combinatorial a leading role in the phase stabilization; in fact, the rp (A) radius approaches. Additionally, formulas from GEN3 show the lowest appears in the best-performing formulas more frequently than the avg(RMSE) among all the others. We also note, from Table I, that other basic atomic properties. Therefore, we further analyzed the GEN1 and GEN3 provide lower avg(RMSE) compared to GEN2 dependence of ΔE on rp (A). In Fig. 2, we show ΔE as a function of and GEN4, respectively; however, GEN2 and GEN4 have a higher rp (A), including fitting curves proportional to rp (A)2 and rp (A)3 . success rate in terms of classification prediction. (For instance, in What can be observed is a clear dependence of ΔE on rp (A): larger Table S2 of the supplementary material, we report the formula (smaller) rp (A) favors RS (ZB). Moreover, there is an overall good optimized using the success rate as a target label.) It is noteworthy agreement with the fit, particularly using the rp (A)3 function. The that the best formula from GEN4 shows similar terms to the corre- latter is, in fact, the most recurrent prototype function detected by sponding case for the ΔE-based optimization. This testifies the fact the ML models. Such a strong dependence of the energy is not that the choice of the performance metrics to rank the material fea- observed with respect to the other atomic properties; other compar- tures can be different according to the target problem to be studied; ative plots of ΔE as a function of other f (p) are reported in Fig. S2 different models’ performance metrics are, in fact, not always of the supplementary material. This behavior is in line with the correlated. further observation that the rocksalt structures systematically show In order to gather hints on the relative contribution of the larger interatomic distances with respect to the zinc blende coun- individual primary atomic properties to the stabilization of either terparts (cf. lattice parameters reported in Ref. 55); therefore, larger J. Appl. Phys. 131, 215703 (2022); doi: 10.1063/5.0088177 131, 215703-7 Published under an exclusive license by AIP Publishing Journal of ARTICLE scitation.org/journal/jap Applied Physics cations prefer to adopt octahedral coordination (i.e., RS) with from the results that the well-performing GEN3 and GEN4 contain longer bond-lengths—and bigger polyhedral volume—compared to all the four radii [rs (A), rs (B), rp (A) and rp (B)], as expected from a ZB with tetrahedral coordination. basic understanding of bonding in octet binary semiconductors. From the obtained results, we remark that formulas based on Note that GEN3 and GEN4 generally show better performance as “spatial” atomic properties achieve higher ranking, thus better per- they are built to explore a wider space of search (see Sec. II and formance, with respect to those, including atomic energy terms, Table S4 in the supplementary material where the number of gen- both in the original models and in the optimized ones. erated and evaluated formulas is also reported). Accordingly, this behavior further confirms the primary role played With the aim of further proving such trends and validate the by the atomic size (in terms of steric and/or bonding-related implemented combinatorial ML method, we study the energetics effects), in determining the energetics of the AB compounds, i.e., in alloys of the type [Ax A01x ]B and A[Bx B01x ], where x is the rela- in selecting the preferred crystal structure.55 In particular, we note tive concentration of the mixing ions, monotonically tuning thus 11 December 2023 10:28:07 FIG. 5. Comparison of actual (i.e., DFT) vs predicted total energy difference ΔE for (a) 1D, (c) 2D, and (d) 3D formulas constructed using GEN3. Panel (b) shows the best 1D descriptors after formula optimization. Lower-right insets show a zoom in the relevant region where many compounds are concentrated. Red dotted lines corre- spond to a 2  avg(RMSE) value. The respective descriptors can be inferred from Tables I–IV. J. Appl. Phys. 131, 215703 (2022); doi: 10.1063/5.0088177 131, 215703-8 Published under an exclusive license by AIP Publishing Journal of ARTICLE scitation.org/journal/jap Applied Physics the average size of one ion with respect to the other. All the alloy We remark that the process of formula optimization is less input properties were linearly interpolated between corresponding computationally expensive than the construction of higher- values for end binaries (i.e., AB and A0 B in the [Ax A01x ]B case) dimensional formulas. In addition, from the formula optimization, according to Vegard’s law.79 For the A-ion mixing case, we con- one can gain better physical insights into the contribution of indi- sidered SrSe, CaSe, MgSe, and BeSe as parent AB compounds, vidual primary atomic properties. These comments overall suggest already included in the original dataset. We then predicted the that lower-dimensional formulas constitute a better choice in terms energy differences between RS and ZB phases for varying concen- of physical interpretation and computational efficiency. trations using the original and optimized 1D formulas constructed via GEN3 and GEN4 generators (Tables I and II, respectively). To IV. CONCLUSIONS confirm the obtained predictions, we thus calculated the energy dif- ference via DFT simulations for a few intermediate concentrations. The knowledge of a material stable crystal structure constitutes The results, shown in Fig. 3, demonstrate an overall agreement the starting point for any ab initio modeling since material proper- between first-principles calculated and machine-learning predicted ties crucially depend on the periodic atomic arrangement in the energetics. In particular, we notice a change of sign in ΔE, reflect- crystal. Within this general framework, our aim here has been to ing the change in the stability of the RS with respect to the ZB exploit ML methods to correlate the energetic stability of different phase, when moving from the larger strontium to the smaller beryl- crystal structures (zinc blende vs rock salt) for popular binary semicon- lium at the A-site, in line with the previously discussed relation ducting compounds with primary properties of their atomic constitu- between atomic radii of the A-ions and phase stabilization. At vari- ents, the latter representing simple and easily accessible ingredients. ance, no such change of phase is observed when mixing ions at the Based on atomic properties, we, therefore, built the material features B-site, keeping fixed the A-type one. This is confirmed by looking at using a combinatorial approach, we trained the machine-learning the energetics in B[Sb1x Px ] and Sr[Se1x Sx ] alloys, shown in Figs. 4 model using the created features over a density-functional-theory (a) and 4(b), respectively. Despite the changing size of the average dataset, and we obtained simple mathematical expressions to quantita- B-site, the two systems preserve the crystal structure adopted by the tively predict the energetic stability of one crystal structure over the parent compounds, i.e., rock salt for the Sr-based compounds and zinc other (i.e., a formula). In addition, we have also introduced an extra blende for the B-based compounds. Such behavior is still in line with step following the linear regression to explore the relative contributions the preferred atomic structure fixed by the ion at the A-site, consis- of individual basic atomic properties. tently with strontium being larger than boron. Qualitative agreement To investigate the performance of the combinatorial approach, we compared our results with a reference paper,55 where the 11 December 2023 10:28:07 between ML-predicted and DFT-calculated energetics is observed again. authors predicted the stability of the crystal structure using an After discussing the results related to 1D models, we now automated feature selection method. We found that our 1D formu- comment about the higher-dimensional formulas. Our best 2D and las constructed using the combinatorial approach achieved a higher 3D formulas from different generators are reported in Tables III accuracy with respect to the reference ones. Furthermore, we also and IV, respectively. learned more about the underlying mechanism from the formula To visualize the performance of the obtained formulas, we optimization, where we found that the stability of RS and ZB reproduce in Fig. 5 the scatterplots of DFT-calculated energies as heavily depends on the rp radius of A-sites. This kind of under- a function of model-predicted energy differences for the best for- standing is, in general, much more difficult to achieve in heavily mulas obtained by GEN3—in terms of avg(RMSE)—for 1D, 1D automated artificial-intelligence methods, such as neural networks, after formula optimization, 2D, and 3D models. From these, one where it is not possible to interpret directly the model results. In can infer the quality of the prediction for the different this respect, our approach based on linear regression allows the approaches: the narrower the area between red lines [representing construction of physical models supported by machine-driven sug- 2  avg(RMSE)], the smaller the error or, equivalently, the more gestions of relevant ingredients; as such, it should be regarded as a reliable the prediction. Notably, this is the case when building methodology offering a huge range of applications in addressing higher dimension formulas. microscopic mechanisms underlying different phenomena, calling In addition, a careful comparison between our results and for extensive investigations in the near future. those reported in the reference paper, Ref. 55, is reported in Table S1 of the supplementary material. In particular, in Fig. S1 of SUPPLEMENTARY MATERIAL the supplementary material, we compared the scatterplot of the 1D See the supplementary material for technical details related to formula from GEN3 and Ref. 55, with bar graphs of errors for indi- LR, DFT calculations of the alloy supercell, dataset, and for addi- vidual compounds. To check the improvement with respect to 1D tional results related to 1D, 2D, and 3D formulas. formulas, we considered the avg(RMSE) value, as also chosen in Ref. 55. One can observe the improvement in avg(RMSE) if we examine 1D and 2D formulas in Tables I and III. We notice ACKNOWLEDGMENTS around 10%–20% improvement from the original 1D to 2D, but This project has received funding from the European Union’s less than 10% of optimized 1D to original 2D formulas. Horizon 2020 research and innovation program under the Marie Furthermore, we also notice that original and optimized 1D formu- Sklodowska-Curie grant (Agreement No. 861145-BeMAGIC). The las from GEN3 and GEN4 better perform with respect to the corre- authors acknowledge the Italian MIUR (Ministry of Education, sponding 2D ones reported in Ref. 55. University and Research) for supporting the PRIN project “TWEET: J. Appl. Phys. 131, 215703 (2022); doi: 10.1063/5.0088177 131, 215703-9 Published under an exclusive license by AIP Publishing Journal of ARTICLE scitation.org/journal/jap Applied Physics 19 ToWards fErroElectricity in Two dimensions” (Grant No. A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, 2017YCTB59) and the “Nanoscience Foundries and Fine Analysis” D. Gunter, D. Skinner, G. Ceder, and K. A. Persson, APL Mater. 1, 011002 (NFFA-MIUR Italy) project. Calculations were performed exploiting (2013). 20 L. Talirz, S. Kumbhar, E. Passaro, A. V. Yakutovich, V. Granata, F. Gargiulo, the computing resources at the Pharmacy Department, Univ. M. Borelli, M. Uhrin, S. P. Huber, S. Zoupanos, C. S. Adorf, C. W. Andersen, Chieti-Pescara. D.A. is grateful to M. Verstraete and B. Dupé (ULiège) O. Schütt, C. A. Pignedoli, D. Passerone, J. VandeVondele, T. C. Schulthess, for the time allowed to work on the writing of this paper. We are also B. Smit, G. Pizzi, and N. Marzari, Sci. Data 7, 299 (2020). thankful to L. Ghiringhelli for his fruitful support and insights. 21 G. Pizzi, A. Cepellotti, R. Sabatini, N. Marzari, and B. Kozinsky, Comput. Mater. Sci. 111, 218 (2016). 22 S. P. Huber, E. Bosoni, M. Bercx, J. Bröder, A. Degomme, V. Dikan, K. Eimre, AUTHOR DECLARATIONS E. Flage-Larsen, A. Garcia, L. Genovese, D. Gresch, C. Johnston, G. Petretto, Conflict of Interest S. Poncé, G.-M. Rignanese, C. J. Sewell, B. Smit, V. Tseplyaev, M. Uhrin, D. Wortmann, A. V. Yakutovich, A. Zadoks, P. Zarabadi-Poor, B. Zhu, The authors have no conflicts to disclose. N. Marzari, and G. Pizzi, npj Comput. Mater. 7, 136 (2021). 23 H. Park, A. Ali, R. Mall, H. Bensmail, S. Sanvito, and F. El-Mellouhi, Mach. Learn.: Sci. Technol. 2, 025030 (2021). DATA AVAILABILITY 24 C. Kim, G. Pilania, and R. Ramprasad, Chem. Mater. 28, 1304 The data that support the findings of this study are available (2016). 25 within the article and its supplementary material. The code for E. Tsymbalov, Z. Shi, M. Dao, S. Suresh, J. Li, and A. Shapeev, npj Comput. Mater. 7, 76 (2021). machine learning is available at https://github.com/lstorchi/ 26 C. J. Bartel, C. Sutton, B. R. Goldsmith, R. Ouyang, C. B. Musgrave, matinformatics. L. M. Ghiringhelli, and M. Scheffler, Sci. Adv. 5, eaav0693 (2019). 27 A. G. Kusne, T. Gao, A. Mehta, L. Ke, M. C. Nguyen, K.-M. Ho, V. Antropov, C.-Z. Wang, M. J. Kramer, C. Long, and I. Takeuchi, Sci. Rep. 4, 6367 REFERENCES (2014). 1 G. E. Moore, Electronics 38, 114–117 (1965). 28 H. Koinuma and I. Takeuchi, Nat. Mater. 3, 429 (2004). 2 29 R. Ramprasad, R. Batra, G. Pilania, A. Mannodi-Kanakkithodi, and C. Kim, npj S. Manti, M. K. Svendsen, N. R. Knøsgaard, P. M. Lyngby, and K. S. Thygesen, Comput. Mater. 3, 54 (2017). “Predicting and machine learning structural instabilities in 2D materials,” 3 M. Fukuda, J. Zhang, Y.-T. Lee, and T. Ozaki, Mater. Adv. 2, 4392 (2021). arXiv:2201.08091 (2022). 4 30 D. Schwalbe-Koda, S. Kwon, C. Paris, E. Bello-Jurado, Z. Jensen, E. Olivetti, K. Kim, L. Ward, J. He, A. Krishna, A. Agrawal, and C. Wolverton, Phys. Rev. 11 December 2023 10:28:07 T. Willhammar, A. Corma, Y. Román-Leshkov, M. Moliner, and Mater. 2, 123801 (2018). R. Gómez-Bombarelli, Science 374, 308 (2021). 31 K. Pal, C. W. Park, Y. Xia, J. Shen, and C. Wolverton, “Scale-invariant 5 E. R. Homer, Comput. Mater. Sci. 161, 244 (2019). machine-learning model accelerates the discovery of quaternary chalcogenides 6 S. Curtarolo, G. L. W. Hart, M. B. Nardelli, N. Mingo, S. Sanvito, and O. Levy, with ultralow lattice thermal conductivity,” arXiv:2109.03751 (2021). Nat. Mater. 12, 191 (2013). 32 M. Kuban, S. Rigamonti, M. Scheidgen, and C. Draxl, “Density-of-states simi- 7 M. L. Green, C. L. Choi, J. R. Hattrick-Simpers, A. M. Joshi, I. Takeuchi, larity descriptor for unsupervised learning from materials data,” arXiv:2201. S. C. Barron, E. Campo, T. Chiang, S. Empedocles, J. M. Gregoire, A. G. Kusne, 02187 (2022). 33 J. Martin, A. Mehta, K. Persson, Z. Trautt, J. V. Duren, and A. Zakutayev, Appl. K. Schütt, P.-J. Kindermans, H. E. Sauceda Felix, S. Chmiela, A. Tkatchenko, Phys. Rev. 4, 011105 (2017). and K.-R. Müller, Adv. Neural Inf. Process. Syst. 30, 992–1002 (2017). 8 A. Walsh, Nat. Chem. 7, 274 (2015). 34 L. Breiman, Mach. Learn. 45, 5 (2001). 35 9 J. Shen, V. I. Hegde, J. He, Y. Xia, and C. Wolverton, Chem. Mater. 33, 9486 K. Gurney, Introduction to Neural Networks (UCL Press Limited, London, (2021). 1997). 10 S. D. Griesemer, L. Ward, and C. Wolverton, Phys. Rev. Mater. 5, 105003 36 T. Xie and J. C. Grossman, Phys. Rev. Lett. 120, 145301 (2018). 37 (2021). R. Guidotti, A. Monreale, S. Ruggieri, F. Turini, F. Giannotti, and 11 J. E. Saal, S. Kirklin, M. Aykol, B. Meredig, and C. Wolverton, JOM 65, 1501 D. Pedreschi, ACM Comput. Surv. 51, 1 (2018). 38 (2013). C. Molnar, Interpretable Machine Learning, 2nd ed. (Christoph Molnar, 2022). 12 S. Kirklin, J. E. Saal, B. Meredig, A. Thompson, J. W. Doak, M. Aykol, S. Rühl, 39 L. Breiman, Mach. Learn. 45, 5 (2001). 40 and C. Wolverton, npj Comput. Mater. 1, 15010 (2015). A. Fisher, C. Rudin, and F. Dominici, arXiv:1801.01489 (2019). 41 13 C. Draxl and M. Scheffler, MRS Bull. 43, 676 (2018). S. Al-Askaar and M. Perkowski, in 2021 IEEE 51st International Symposium 14 C. Draxl and M. Scheffler, J. Phys. Mater. 2, 036001 (2019). on Multiple-Valued Logic (ISMVL) (IEEE, 2021), pp. 128–135. 15 S. Curtarolo, W. Setyawan, S. Wang, J. Xue, K. Yang, R. H. Taylor, 42 R. Elshawi, M. H. Al-Mallah, and S. Sakr, BMC Med. Inform. Decis. Mak. 19, L. J. Nelson, G. L. Hart, S. Sanvito, M. Buongiorno-Nardelli, N. Mingo, and 1 (2019). O. Levy, Comput. Mater. Sci. 58, 227 (2012). 43 M. T. Ribeiro, S. Singh, and C. Guestrin, “Model-agnostic interpretability of 16 M. N. Gjerding, A. Taghizadeh, A. Rasmussen, S. Ali, F. Bertoldo, machine learning,” arXiv:1606.05386 (2016). T. Deilmann, N. R. Knøsgaard, M. Kruse, A. H. Larsen, S. Manti, T. G. Pedersen, 44 S. Lundberg and S.-I. Lee, “A unified approach to interpreting model predic- U. Petralanda, T. Skovhus, M. K. Svendsen, J. J. Mortensen, T. Olsen, and tions,” arXiv:1705.07874 (2017). 45 K. S. Thygesen, 2D Mater. 8, 044002 (2021). Y. Zhang, P. Tiňo, A. Leonardis, and K. Tang, IEEE Trans. Emerg. Top. 17 S. Haastrup, M. Strange, M. Pandey, T. Deilmann, P. S. Schmidt, Comput. Intell. 5, 726–745 (2021). 46 N. F. Hinsche, M. N. Gjerding, D. Torelli, P. M. Larsen, A. C. Riis-Jensen, Q. Zhang, Y. N. Wu, and S.-C. Zhu, in Proceedings of the IEEE Conference on J. Gath, K. W. Jacobsen, J. Jørgen Mortensen, T. Olsen, and K. S. Thygesen, 2D Computer Vision and Pattern Recognition (IEEE, 2018), pp. 8827–8836. 47 Mater. 5, 042002 (2018). D. Alvarez Melis and T. Jaakkola, arXiv:1806.07538 (2018). 48 18 F. Bertoldo, S. Ali, S. Manti, and K. S. Thygesen, “Quantum point defects in S. Chatterjee and J. S. Simonoff, Handbook of Regression Analysis (Wiley, 2D materials: The QPOD database,” arXiv:2110.01961 (2021). Hoboken, NJ, 2013). J. Appl. Phys. 131, 215703 (2022); doi: 10.1063/5.0088177 131, 215703-10 Published under an exclusive license by AIP Publishing Journal of ARTICLE scitation.org/journal/jap Applied Physics 49 63 J. Leskovec, A. Rajaraman, and J. D. Ullman, Mining of Massive Datasets S. Shalev-Shwartz and S. Ben-David, Understanding Machine Learning: From (Stanford University, Stanford, CA, 2010). Theory to Algorithms (Cambridge University Press, New York, 2014). 50 T. Miller, Artif. Intell. 267, 1 (2019). 64 G. Van Rossum and F. L. Drake, Jr., Python Tutorial (Centrum voor 51 B. Kim, R. Khanna, and O. O. Koyejo, Adv. Neural Inf. Process. Syst. 29, Wiskunde en Informatica, Amsterdam, 1995). 2280–2288 (2016). 65 L. Storchi, “Open Source Code”; see https://github.com/lstorchi/matinformatics 52 D. G. Kleinbaum, K. Dietz, M. Gail, M. Klein, and M. Klein, Logistic (2022). 66 Regression (Springer, 2002). MATHEMATICA, version 13.0.0, W. R., Inc., Champaign, IL, 2021. 53 A. J. Myles, R. N. Feudale, Y. Liu, N. A. Woody, and S. D. Brown, J. Chemom. 67 L. Storchi, “Mathematica Notebook,” startingtest.nb; see https://github.com/ 18, 275 (2004). lstorchi/matinformatics (2022). 54 L. M. Ghiringhelli, “Interpretability of machine-learning models in physical 68 F. Gao and L. Han, Comput. Optim. Appl. 51, 259 (2012). 69 sciences,” arXiv:2104.10443 (2021). G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed., Johns 55 L. M. Ghiringhelli, J. Vybiral, S. V. Levchenko, C. Draxl, and M. Scheffler, Hopkins Studies in the Mathematical Sciences (The Johns Hopkins University Phys. Rev. Lett. 114, 105503 (2015). Press, Baltimore, MD, 2013). 56 B. Meredig, A. Agrawal, S. Kirklin, J. E. Saal, J. W. Doak, A. Thompson, 70 C. G. Broyden, IMA J. Appl. Math. 6, 76 (1970). K. Zhang, A. Choudhary, and C. Wolverton, Phys. Rev. B 89, 094104 (2014). 71 R. dembo, S. Eisenstat, and T. Steihaug, SIAM J. Numer. Anal. 19, 400 (1982). 57 L. Ward, A. Dunn, A. Faghaninia, N. E. R. Zimmermann, S. Bajaj, Q. Wang, 72 G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993). J. Montoya, J. Chen, K. Bystrom, M. Dylla, K. Chard, M. Asta, K. A. Persson, 73 G. Kresse and J. Furthmüller, Comput. Mater. Sci. 6, 15 (1996). G. J. Snyder, I. Foster, and A. Jain, Comput. Mater. Sci. 152, 60 (2018). 74 G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996). 58 L. Ward, A. Agrawal, A. Choudhary, and C. Wolverton, npj Comput. Mater. 2, 75 V. Blum, R. Gehrke, F. Hanke, P. Havu, V. Havu, X. Ren, K. Reuter, and 16028 (2016). M. Scheffler, Comput. Phys. Commun. 180, 2175 (2009). 59 P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964). 76 C. Eckhardt, K. Hummer, and G. Kresse, Phys. Rev. B 89, 165201 (2014). 60 W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965). 77 L. Bellaiche and D. Vanderbilt, Phys. Rev. B 61, 7877 (2000). 61 J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992). 78 D. Amoroso, A. Cano, and P. Ghosez, Phys. Rev. B 97, 174108 (2018). 62 K. Momma and F. Izumi, J. Appl. Crystallogr. 41, 653 (2008). 79 L. Vegard, Z. Phys. 5, 17 (1920). 11 December 2023 10:28:07 J. Appl. Phys. 131, 215703 (2022); doi: 10.1063/5.0088177 131, 215703-11 Published under an exclusive license by AIP Publishing