Towards machine learning for microscopic mechanisms: a formula search for crystal structure stability based on atomic properties

Machine Learning (ML) techniques are revolutionizing the way to perform efficient materials modeling. Nevertheless, not all the ML approaches allow for the understanding of microscopic mechanisms at play in different phenomena. To address the latter aspect, we 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 semiconducting binary compounds with respect to zincblende and rocksalt crystal structures. The adopted models are trained using 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 performances at a very small computational cost. In addition, our approach allows 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 stabilization of one crystal structure with respect to the other, suggesting the major relevance of the radius associated to the $p$-shell of the cation species.


I. INTRODUCTION
Modeling material properties with high accuracy and low computational cost is one of the grand-challenges in materials science and engineering. The development of ab-initio methods have provided accurate tools for material properties prediction and their further optimization; nevertheless, one disadvantage of approaches relying only on first-principles simulations is the high cost required in terms of computational resources and simulation time. In recent years, the continuous growth of available computational power 1 has stimulated scientists to move in the direction of high-throughput simulations [2][3][4][5][6][7][8][9][10] . Along this line, open access databases, such as OQMD 1112 , NOMAD 13,14 , 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 basic first-principles results. A large amount of ab-initio data is thus available, which can be used for deeper analyses and studies, provided one can count on proper tools to extract relevant information out of them. Therefore, in the last years, materials scientists have developed different Machine Learning (ML) methods to rationalize the data analysis [23][24][25][26][27][28][29][30][31][32] . Each method has its own specific advantages and limitations. Methods like Neural Network (NN) 33 or Random Forest 34 are very efficient 35 but not always transparent, blurring the comprehension of the role played by the input variables in the final results; ML methods, based for instance on linear regression (LR) 36,37 , appear to be more suitable to obtain predictive and comprehensible models 38,39 . Nevertheless, finding a linear dependence between input and output properties is not always an easy task.
In this work, we thus propose a ML-based approach to build sets of features (or descriptors) starting from a given set of basic variables (e.g., atomic properties), which are subse-quently used to construct LR models (or formulas). To test our method, we target a prototypical case in material science: the classification of the most stable crystal structure between rock-salt (RS) and zinc-blende (ZB) for semiconductor AB binary compounds 40 . In our approach, we adopt both simple one-dimensional and multi-dimension LR. To identify useful features, we generate combinations of basic atomic properties (i.e. the independent variables in our approach) of the material constituents through a combinatorial approach 41 . We then carry out an analysis of the emerging best-performing formulas, identifying the role of specific atomic features in determining the final stabilization of the crystal structure. Finally, we test the predictive capability of the obtained formulas by applying them to "new" compounds (i.e. outside the dataset used for training the model), finding an overall satisfactory agreement with first-principles results. We remark that our approach is similar to what originally proposed by Ghiringhelli et al. 40 , though with some differences and further extensions, which will be carefully discussed in what follows.

II. METHODOLOGY
The approach we present here can be regarded as a combinatorial machine-learning: a set of basic atomic properties (APs, listed in Table S.2 in Supplementary Information) are randomly combined (though under certain initial constraints detailed below), to build a set of material features (MFs). The generated features are then used to train a LR model, where the energy difference between rocksalt and zincblende structures is the dependent variable (i.e., the label). Then, we select the best performing model according to standard performances metrics, such as the Root Mean Squared Error (RMSE). The final result of this procedure is a "formula", which is a concise and clear representation of the relationship between the used atomic properties and the energy difference between RS and ZB phases. In the following, we describe in detail the different steps of our approach.

A. Dataset preparation and materials
As mentioned, we aim at predicting the total energy difference (∆E = E RS − E ZB ) between RS and ZB phases of cubic crystal structures for 82 semiconductor binary AB compounds (the dataset is reported in table S.2 in SI). We employed total energies reported in Ref. 40 , which were calculated through density functional theory (DFT) 4243 within the local density approximation (LDA 44 ).
The construction of the material features (MFs), is based on primary atomic properties of the constituents, also taken from Ref. 40 . To facilitate the physical interpretation of each MF, the APs are subdivided into two different kinds: (i) "energy" properties, including highest occupied Kohn-Sham level (HOMO), lowest unoccupied Kohn-Sham level (LUMO), Ionization Potential (IP), Electron affinity (EA); (ii) "spatial" properties, including r s , r p , and r d , i.e. the radii where the radial probability density of the valence s, p, and d orbitals, respectively, reaches its maximum.

Formulas construction
We rely on the LR 3637 approach to obtain a direct interpretation of the dependent and independent variables. The construction of a useful LR model can become troublesome, requiring a linear dependence between features. In Ref. 40 , the authors implemented an automated feature selection method employing the LASSO regression analysis method 40,46 . In our work, we use a combinatorial approach to generate the dependent variable (material features) to be used within the linear equations, and thus to finally obtain the formulas.
In Fig. 1, we illustrate the workflow of the formula generation and selection using LR. The process starts with the selection of the APs to be combined. Afterwards, we choose prototype functions that are simple analytical operations applied to the APs. In our case, we selected 5 prototype functions, f (x), namely x, x 2 , x 3 , √ x, e x . where x is an AP. Then, we obtain the final set of MFs by combining different prototype functions via the combinatorial approach (see for instance 41 ), and applying the following additional set of rules: • GEN1: combine two prototype functions in the numerator, forcing them to belong to the same kind of APs, that is both "spatial"-like or both "energy"-like; one prototype function is at the denominator with the only constraint to be non-zero, such as • GEN2: combine two prototype functions with same kind of APs at the numerator, and a single prototype function at the denominator with argument of a different kind with respect to the numerator ones. For instance, if AP 1 in f 1 (AP 1 ) and AP 2 in f 2 (AP 2 ) is an "energy" term (i.e. EA or HOMO), then AP 3 must be a "spatial" term, (i.e. r p ) • GEN3: combine two prototype functions at both the numerator and denominator without any constraints • GEN4: combine two prototype functions with the same physical dimensions at both the numerator and denominator where = + − ×÷ Each one of these set of rules corresponds to a different MFs generator.
From the implementation point of view, each generator is a Python 47 function that produces a set of strings. Therefore, we can easily exploit the Python capability to parse a source code and run Python expression (code) within a program 48 to compute all the MFs' values starting from the generated sets of strings. This allows for an easy implementation and plugin of other generators, leaving the workflow unchanged: a new generator can be introduced implementing a Python function returning a list of strings, each one being a valid MF.
Finally, in order to choose the optimal formula, we build a LR model for each of the generated MF. To practically select the best model, i.e. the "best formula", we randomly split the full dataset into : 90% as training set to train/initialize the model; 10% as a test set to check model's performance. We perform this random splitting N times (with N = 150) for each model, and we calculate the RMSE from the test set for each run. Afterwards, we again verify the top 10 resulting best formulas with a higher value of training set and test set splitting, with N = 1000. We average it over all N splitting, and we obtain avg(RMSE), as reported in our Tables.
We mention that different metrics for evaluating regression model can lead to different formulas ranking. In this work, we rank the obtained models based on the lowest avg(RMSE) for direct comparison with a previous work 40 .

B. Formula optimization
In order to further improve the performance of our models, we introduce an additional step, which we refer to as "formula optimization". In detail, we focus on the top 10 formulas obtained using each generator and the subsequent LR, as described in the previous section. After that, we use a grid  45 . Grey (yellow) spheres represent A (B) atoms. c) Workflow for formulas construction, machine learning methodology, validation, and MF selection procedures. In the AB compounds, A is the atom with the lowest electronegativity.
search to find the relative weights of each prototype function of the atomic properties (i.e., each f i (AP i )) within the formula. A first grid search ranging between -1 to 1 with the increasing step of 0.1 is used. We multiply each f i (AP i ) of the formula by the weight coefficient and we optimize the final RMSE value. Once the procedure finds a set of optimal weight coefficients, two subsequent grid searches, with reduced incremental step values (0.01 and 0.001 respectively) and range of search are performed to obtain the final set of refined weight coefficients. Noteworthy, for each set of weight coefficients generated during the grid search, we also run the linear regression. Thus, we are performing a proper formula optimization, as at each step of the grid search we are updating both the weight coefficients as well as the slope and intercept coming from the LR.
To further clarify the procedure, we show here an exemplary equation: where ∆E is the targeted material feature (MF), a, b, c, d denote the weight coefficients scanned during the grid search, are the prototype functions build on the primary atomic properties AP i , and m and q are the the slope or angular coefficient and intercept, respectively, recursively determined upon LR.
In Table II, we report the optimized, best performing formula from the different generators; the top 10 formulas are reported in Table S.1 of the Supplementary Material.
To benchmark our grid search, we also used automated coefficient-optimizing methods: Nelder-Mead 49 , Conjugate Gradient (CG) 50 , Broyden-Fletcher-Goldfarb-Shanno (BFGS) 51 and Truncated Newton method (TNC) 52 . Although the resulting sets of coefficients are different in terms of single values with respect to those obtained via the grid search, the ratios between them is almost preserved as well as the associated RMSE. In particular, for the case of GEN1 and GEN2, the ratio between the numerator coefficients a and b is preserved; for GEN3 and GEN4 also the denominator coefficients ratio, between c and d, is preserved. In Fig. S.3 of the Supplementary Material, we show the evolution of the RMSE and different ratios for different methods using 1D feature generated by GEN3.

C. Higher-dimensional features
For the construction of higher dimensional 2D-formulas, we combined in all possible ways pairs of MFs extracted from the best 1000 ones and checked the avg(RMSE) using multiple LR for N test-train set splits. We followed the same process to construct the 3D formulas, where three different 1D MFs are combined. The comparison between performances is discussed in the following Section.

D. Test of predictive power of ∆E formula for novel AB compounds
After obtaining the optimised 1D formulas for ∆E in the case of AB compounds, we aimed at further verifying their validity and predictive power, by considering additional AB systems (i.e. which were not originally included in the ML training set) and by comparing values obtained from ML-predicted ∆E formula with corresponding ab-initio calculated values. In closer detail, we focused on different alloys, obtained by changing respectively the concentration of A-site atoms, such as [A x A 1−x ]B, and of B-site atoms, such as A[B x B 1−x ]. Accordingly, one can test the efficiency of the formulas by check-ing the energy difference for intermediate concentrations as obtained from optimised 1D formulas and compare their trend with respect to first-principles results. To this end, ab-initio electronic-structure simulations were carried out within DFT and LDA functional. Calculations were performed using the VASP 53-55 code, employing a 8 × 8 × 8 k-mesh for the Brillouin zone sampling. We verified that the results obtained with the pseudopotential VASP for the parent binary compounds were consistent with those reported by Ghiringhelli et al., calculated with the all-electron FHI-aims code 56 . For simulations at different concentrations, we adopted the so called "Virtual Crystal Approximation" (VCA), based on virtual atoms interpolating between the real constituent atoms 57,58 . However, as well known from the literature, 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 alloys properties 59 . Accordingly, in some cases (i.e. for Mg x Ca 1−x Se alloys), in order to mimic disordered structures with an improved accuracy, we calculated total energies using supercell structures, rather than using the VCA method on primitive unit cells. Specifically, the considered supercell is the cubic unit cell composed by four AB formula units with planes of cations alternating along the c direction (see Figure S.4). The k-mesh was modified accordingly, to maintain the same density of points employed in the simulations of primitive cells.

III. RESULTS AND DISCUSSION
In this section, we will analyse the final formulas as obtained from different generators. The results are shown in Tables I,II,III,IV; in the first row we report the results obtained by Ghiringhelli et al. 40 for comparison.
First, by comparing the avg(RMSE) values, we note that all 1D formulas obtained from our different generators better perform with respect to the 1D ones reported in 40 , where the authors used the automated feature selection method LASSO 46 . Noteworthy, some atomic primary features appearing in 1D formulas of Ref. 40 also appear in our obtained list of 1D formulas using GEN1 and GEN2; nevertheless, those are characterized by a higher avg(RMSE) than other formulas we obtained via our combinatorial approaches. Additionally, formulas from GEN3 show the lowest avg(RMSE) among all the others. We also note, from Table I, that GEN1 and GEN3 provide lower avg(RMSE) compared to GEN2 and GEN4 respectively; however, GEN2 and GEN4 have a higher success rate in terms of classification prediction. This testifies the fact that the choice of the performances metric to rank the material features can be different according to the target problem to be studied; different models' performances metric are, in fact, not always correlated.
In order to gather hints on the relative contribution of the individual primary atomic properties to the stabilization of either the rocksalt or the zincblende structure, we extracted the best ten formulas with the lowest avg(RMSE) from each generator (so called "original" formulas) and then apply the formula optimization, as detailed in the previous section. This procedure attributes relative weights to each f (AP), allowing to measure the importance of the individual atomic properties in driving the energy stabilization. In principle, the avg(RMSE) value depends on random test-train splits that we perform to our dataset. Therefore, to reduce the effect of randomization, as a target model performances metric, we rank our optimized formula based on the RMSE of the whole dataset, rather than based on avg(RMSE). By comparing Table I and Table II, it is evident that the optimization procedure can further change the formulas ranking, providing a different final "best formula" with respect to the non-optimized formulas. In particular, we notice an improvement in RMSE around 5-10% after the formula optimization.
Interestingly, our results reveal the size of the A-ion to play a leading role in the phase stabilization; in fact, the r p (A) radius appears in the best performing formulas more frequently than the other basic atomic properties. Therefore, we further analysed the dependence of ∆E on r p (A). In Fig. 2,  From the obtained results, we remark that formulas based on "spatial" atomic properties achieve higher ranking, thus better performance, with respect to those including atomic energy terms, both in the original models and in the optimized ones. Accordingly, this behaviour further confirms the primary role played by the atomic size (or, equivalently, steric effects), in determining the energetics of the AB compounds, i.e. in selecting the preferred crystal structure 40 . In the aim of further proving such trends and validate the implemented combinatorial ML method, we study the energetics in alloys of the type , where x is the relative concentration of the mixing ions, monotonically tuning thus the average size of one ion with respect the other. All the alloy input properties were linearly interpolated between corresponding values for end binaries (i.e. AB and A B in the [A x A 1−x ]B case), according to the Vegard's law 60 . For the A-ion mixing case, we considered SrSe, CaSe, MgSe, BeSe as parent AB compounds, already included in the original dataset. We then predicted the energy differences between RS and ZB phases for varying concentrations using the original and optimized 1D formulas constructed via GEN3 and GEN4 generators (Table I and  Table II, respectively). To confirm the obtained predictions, we thus calculated the energy difference via DFT simulations, for a few intermediate concentrations. The results, shown in Fig. 3, demonstrate an overall agreement between first-principles calculated and machine-learning predicted energetics. In particular, we notice a change of sign in ∆E, reflecting the change in the stability of the RS with respect to the ZB phase, when moving from the larger Strontium to the smaller Beryllium at the A-site, in line with the previously discussed relation between atomic radii of the A-ions and phase stabilization. At variance, no such change of phase is observed when mixing ions at the B-site, keeping fixed the A-type one. This is confirmed, by looking at the energetics in B[Sb 1−x P x ] and Sr[Se 1−x S x ] alloys, shown in Fig.4(a) and Fig. 4(b), respectively. Despite the changing size of the average B-site, the two systems preserve the crystal structure adopted by the the parent compounds, i.e. rocksalt for the Sr-based compounds and zincblende for the B-based compounds. Such a behavior is still in line with preferred atomic structure fixed by the ion at the A-site, consistently with Strontium being larger than Boron. Qualitative agreement between ML-predicted and DFT-calculated energetics is observed again.
After discussing the results related to 1D models, we now comment about the higher dimensional formulas. Our best 2D and 3D formulas from different generators are reported in Tables III and IV, respectively. To visualize the performance of the obtained formulas, we reproduce in Fig. 5 the scatter plots of DFT-calculated energies as a function of model-predicted energy differences for the best formulas obtained by GEN3 -in terms of avg(RMSE) -for 1D, 1D after formula optimization, 2D, and 3D models. From these, one can infer the quality of the prediction for the different approaches: the narrower the area between red lines (representing 2×avg(RMSE)), the smaller the error or, equivalently, the more reliable the prediction. Notably, this is the case when building higher dimension formulas.
In addition, a careful comparison between our results and those reported in the reference paper, Ref. 40 , is reported in Table S.1 of the Supplementary Material. In particular, in Fig. S.1 we compared the scatter plot of the 1D formula from GEN3 and Ref. 40 , with bar graphs of errors for individual compounds. To check the improvement with respect to 1D formulas, we considered the avg(RMSE) value, as also chosen in Ref. 40 . One can observe the improvement in avg(RMSE) if we examine 1D and 2D formulas in Tables I and III. We notice around 10-20% improvement from the original 1D to 2D, but less than 10% of optimized 1D to original 2D formulas. Furthermore, we also notice that original and optimized 1D formulas from GEN3 and GEN4 better perform with respect the corresponding 2D ones reported in Ref 40 .
We remark that the process of formula optimization is less computationally expensive than the construction of higherdimensional formulas. In addition, from the formula optimization one can gain a better physical insights about the contribution of individual primary atomic properties. These comments overall suggest that lower-dimensional formulas constitute a better choice in terms of physical interpretation and computational efficiency.

IV. CONCLUSIONS
The knowledge of a material stable crystal structure constitutes the starting point for any ab-initio modelling, since materials properties crucially depend on the periodic atomic arrangement in the crystal. Within this general framework, our aim here has been to exploit ML methods to correlate the energetic stability of different crystal structures (zincblende vs rocksalt) for popular binary semiconducting compounds with primary properties of their atomic constituents, the latter representing simple and easily-accessible ingredients. Based on atomic properties, we therefore built the material features using a combinatorial approach, we trained the machine learning model using the created features over a density-functionaltheory dataset and we obtained simple mathematical expressions to quantitatively predict the energetic stability of one crystal structure over the other (i.e., a formula). In addition, we have also introduced an extra step following the linear regression to explore the relative contributions of individual basic atomic properties.
To investigate the performance of the combinatorial approach, we compared our results with a reference paper 40 , where the authors predicted the stability of the crystal structure using an automated feature selection method. We found that our 1D formulas constructed using the combinatorial approach achieved a higher accuracy with respect to the reference ones. Furthermore, we also learned more about the underlying mechanism from the formula optimization, where we found that the stability of RS and ZB heavily depends on the r p radius of A-sites. This kind of understanding is, in general, much more difficult to achieve in heavily-automated artificialintelligence methods, such as neural networks, where it is not possible to interpret directly the model results. In this respect, our approach based on linear regression allows the construction of physical models supported by machine-driven suggestions of relevant ingredients; as such, it should be regarded as a methodology offering a huge range of applications in addressing microscopic mechanisms underlying different phenomena, calling for extensive investigations in the nearby future.
Towards machine learning for microscopic mechanisms: a formula search for crystal structure stability based on atomic properties 6 V. AUTHOR DECLARATIONS

A. Conflict of Interest
The authors have no conflicts to disclose.

B. Data Availability
The data that supports the findings of this study are available within the article and its supplementary material.