Automated machine learning structure-composition-property relationships of perovskite materials for energy conversion and storage

Dataset of perovskite materials

The dataset of ABO₃ perovskite oxide utilized in this paper originates from DFT high-throughput compute, including the ion radius of the A-site, the ion radius of the B-sit, formation energy, crystal volume, band gap, lattice parameters (a, b, c, α, β, γ), oxygen vacancy formation energy, stability, etc.^[30]. Different radii of ions were used in A sites and B sites, including those of A ∈ [Al, As, Ag, Be, B, Bi, Ba, Ca, Co, Cu, Cr, Cd, Ce, Zr, Zn, Dy, Er, Fe, Ge, Gd, Ga, Hf, Ho, In, Ir, K, La, Lu, Mn, Mg, Mo, Ni, Nb, Nd, Pr, Pd, Ru, Pb, Rb, Re, Rh, Si, Sc, Sr, Sb, Sn, Sm, Ta, Th, Tb, Te, Ti, Tm, U, V, W, Y, Yb, etc.] and B ∈ [Al, Ag, As, Li, B, Bi, Be, Ca, Zn, Co, Cu, Cd, Cr, Ce, Zr, Dy, Eu, Er, Fe, Ga, Ge, Gd, Hf, Ho, In, Ir, Lu, Mn, Mg, Mo, Ni, Nd, Nb, Pr, Pb, Pd, Ru, Rb, Rh, Re, Si, Sb, Sc, Sn, Sm, Sr, Ti, Ta, Th, Tb, Te, Tm, U, V, W, Y, Yb, etc.], for perovskite ABO₃ compounds. Through a cursory examination of the data, a total of 4912 sets of ABO₃ perovskite compound high-throughput data were chosen. The values of stability were in the -0.729~3.927 eV/atom range. Data sets are indexed via abbreviations to make experimentation easier, using the following details: r_A for ionic radius at the A-site, r_B for ionic radius at the B-site, ΔH_f for formation energy, V for crystal volume, ΔE for band gap, a for lattice parameter a, b for lattice parameter b, c for lattice parameter c, α for the lattice parameter of the α phase, β for the lattice parameter of the β phase, γ for the lattice parameter of the γ phase, for oxygen vacancy formation energy and ΔH_s for stability. For ABO₃ perovskite oxide, more detailed descriptions are shown in Table 1.

Table 1

Explanation for ABO3 perovskite compounds datasheet

No.	Property	Unite	Description
1	rA	Å	Ionic radius of A site
2	rB	Å	Ionic radius of B site
3	ΔHf	eV/atom	Formation energy as calculated by equation of the distortion with the lowest energy
4	ΔHs	eV/atom	Stability as calculated by equation of the distortion with the lowest energy
5	V	Å3/atom	Volume per atom of the relaxed structure
6	ΔE	eV	PBE band gap obtained from the relaxed structure
7	a	Å	Lattice parameter a of the perovskite structure
8	b	Å	Lattice parameter b of the perovskite structure
9	c	Å	Lattice parameter c of the perovskite structure
10	α		α angle of the relaxed structure. α = 90 for the cubic, tetragonal and orthorhombic distortion
11	β		β angle of the relaxed structure. β = 90 for the cubic, tetragonal and orthorhombic distortion
12	γ		γ angle of the relaxed structure. γ = 90 for the cubic, tetragonal and orthorhombic distortion
13		eV per O atom	Oxygen vacancy formation energy

Data pretreatment

Data pretreatment is used to process the missing and repeated values in the data, raising the data’s accuracy and helping to raise the precision and efficiency of the subsequent learning procedure. The common processes of data pretreatment include missing value processing, attribute coding, feature selection, etc.^[31,32]. There are three common ways to deal with missing values: use the feature that contains the missing value directly, delete the feature that contains the missing value (this only works if the feature contains blank values in a big number), and complete the missing value^[33,34]. Because there are a small number of blank values in the raw Dataset, the features containing blank values are employed in this paper to process blank values.

Feature selection refers to the procedure of picking a subset of relevant features from a given feature collection^[35]. Although a variety of factors influence the target characteristics of perovskite oxides, the amount of features must be appropriate, the features must be uneven for the category of interest, and certain non-essential information must be removed^[36]. Correlation is a term that describes the degree and direction of the link between these two measurable features. Pearson correlation analysis is often adopted for analyzing the connections between two measurable features^[37]. In this paper, Pearson correlation coefficient was used to examine the link between composition, structure and property, and the linear correlation among composition, structure and property is measured^[38]. Pearson correlation coefficient can be defined easily as follows:

where x_i represents the value of feature x, y_i represents the value of feature y, represents the average value of value of feature x, represents the average value of feature y, σ_x represents the standard deviation of feature x, σ_y represents the standard deviation of feature y and n represents the sample size^[38]. Correlation coefficient of 1 indicates strongly positively correlated, whereas a correlation value of -1 indicates a strong negative correlation. The correlation coefficient near to 0 implies that there is no association^[37,38].

Pearson correlation coefficients were utilized to choose the raw dataset in this paper. Table 2 depicts the degree of connection between the 12 properties of perovskite oxide ABO₃ and their stability. Figure 3 depicts Pearson correlation map for different features.

Table 2

Pearson correlation coefficients of built-up features and stability

PCC	rA	rB	ΔHf	V	ΔE	a	b	c	α	β	γ		ΔHs
rA	1.000000	0.008855	-0.364382	0.183594	0.088459	0.286041	0.272757	0.257065	0.012843	0.012120	0.013025	0.371450	-0.522209
rB	0.008855	1.000000	0.023527	0.754276	0.055515	0.384344	0.370947	0.295416	-0.045394	-0.047375	-0.045085	-0.517963	0.098718
ΔHf	-0.364382	0.023527	1.000000	-0.140423	-0.472760	-0.448701	-0.450275	-0.431603	0.095894	0.094795	0.095947	-0.424975	0.594845
V	0.183594	0.754276	-0.140423	1.000000	0.008301	0.317787	0.277428	0.150485	-0.012753	-0.012908	-0.012622	-0.242627	0.127233
ΔE	0.088459	0.055515	-0.472760	0.008301	1.000000	0.337769	0.332856	0.340382	-0.126428	-0.124734	-0.126561	0.169527	-0.332567
a	0.286041	0.384344	-0.448701	0.317787	0.337769	1.000000	0.982619	0.879740	-0.366276	-0.365131	-0.366207	-0.076172	-0.591562
b	0.272757	0.370947	-0.450275	0.277428	0.332856	0.982619	1.000000	0.913307	-0.305305	-0.305172	-0.305283	-0.079756	-0.610395
c	0.257065	0.295416	-0.431603	0.150485	0.340382	0.879740	0.913307	1.000000	-0.054803	-0.055040	-0.054817	-0.065271	-0.635102
α	0.012843	-0.045394	0.095894	-0.012753	-0.126428	-0.366276	-0.305305	-0.054803	1.000000	0.996639	0.999847	0.062728	0.208962
β	0.012120	-0.047375	0.094795	-0.012908	-0.124734	-0.365131	-0.305172	-0.055040	0.996639	1.000000	0.995354	0.062613	0.207678
γ	0.013025	-0.045085	0.095947	-0.012622	-0.126561	-0.366207	-0.305283	-0.054817	0.999847	0.995354	1.000000	0.062694	0.209099
	0.371450	-0.517963	-0.424975	-0.242627	0.169527	-0.076172	0.079756	-0.065271	0.062728	0.062613	0.062694	1.000000	-0.359896
ΔHs	-0.522209	-0.522209	0.594845	0.127233	-0.332567	-0.591562	-0.610395	-0.635102	0.208962	0.207678	0.209099	-0.359896	1.000000

Figure 3. This is a Pearson correlation map for raw data. The correlation coefficient is shown by the color bar: red indicates strongly positive correlations, white denotes strongly negative correlations. The worth of the related Pearson correlation coefficient is represented by the filled fraction in each tiny square.

Extreme feature engineering

For the sake of rapidly discovering the connections between structure, composition, and properties, we displayed dataset’s feature distribution. Figure 4 depicts the distribution of raw features and stability. The distribution of the raw data set of observation feature and the predicted variable stability is positively biased, and the range of data is quite broad. Therefore, data transformation methods must be employed to generate new descriptors through feature engineering^[39]. Essentially, the data provided to the algorithm should be compatible with the required structure or characteristics of the underlying data. Feature engineering is the process of turning data attributes into data features and extracting features from raw data through algorithms and models to the greatest extent possible^[40]. Therefore, the feature engineering approach may generate a large number of new descriptors and assess their performance with a subset of them.

Figure 4. (A) The radius distribution of A-site ions; (B) the radius distribution of B-site ions; (C) distribution of formation energy; (D) distribution of band gap feature; (E) distribution of volume; (F) distribution of lattice constant a; (G) distribution of lattice constant b; (H) distribution of lattice constant c; (I) distribution of α angle of the crystal structure; (J) distribution of β angle of the crystal structure; (K) distribution of γ angle of the crystal structure; (L) distribution of stability; and (M) distribution of oxygen vacancy formation energy.

In machine learning, feature engineering is a critical data preparation activity that creates suitable descriptors from a given feature to improve prediction performance^[41]. Feature engineering is adding some functions of conversion, such as arithmetic and aggregation operators, into a given attribute to create a huge number of new descriptors^[42]. The transformation functions contribute to increase the dimensions of features or to turn the nonlinear connection between features and stability into a more understandable linear one^[40,43]. Feature combination is a highly important method in feature engineering to integrate features from several categories into a single feature^[44]. This is a beneficial method when a combination of features outperforms a single feature. The feature combination is the cross multiplication of all conceivable eigenvalues in mathematics. The features of each combination really constitute the information synergy.

A huge number of brand-new descriptors were gained through extreme feature engineering, where the dimensionality of features was also expanded. Figure 5 shows the construction process of the descriptors by extreme feature engineering. In their midst, x_i (i = 1, 2, … n) indicates the selected feature. The parameters following the yellow arrows reflect a significant number of new descriptors that were produced. These 9 functions of x, x^-1, , x², x³, e^x, ln|x|, ln(1 + |x|) and log|x| are utilized for nonlinear transformation of features. In order to generate additional descriptors, these descriptors are merged non-linearly^[45]. The primary descriptors were generated in the following way:

Figure 5. The process of descriptor construction by extreme feature engineering. The x_i (i = 1, 2, … n) indicates the chosen feature. The parameters following the arrows denote the constructed a large number of brand-new descriptors. These 9 functions of x, x^-1, , x², x³, e^x, ln|x|, ln(1 + |x|) and log|x| are employed in the nonlinear conversion of features. All of the descriptors mixed in a nonlinear way to construct more descriptors.

• Step 1: Import the chosen vital features into these 9 functions of x, x^-1, , x², x³, e^x, ln|x|, ln(1 + |x|) and log|x|, where x is one of the vital features chosen from the raw features, and it can directly generate brand-new descriptors;

• Step 2: Feature combination combines the brand-new descriptors from Step 1. Increase the number of descriptive words by multiplying them by two or more and then combining them into a brand-new one;

• Step 3: Substituting brand-new descriptors gained in step 2 into the function x^-1 for nonlinear conversion, and the number of descriptors acquired has increased.

Regression

The term “regression” refers to the process of determining the quantitative connection between two or more variables using a group of data, the establishment of simulations from mathematics, and the estimation of unidentified factors^[46]. Machine learning is an efficient way of performing regression. The capacity to do linear regression is defined as to properly depict the connection between data using a straight line, which is more suited to fitting the expression^[47]. The modeling speed of linear regression is rapid, it does not need sophisticated calculation, and it may even run quickly when dealing with huge amounts of data^[48]. The gained linear expression can be understood and interpreted according to the coefficient of each variable, and the influence of each feature on the result can be directly seen from the weight, which is much easier to grasp^[43,49]. Nonlinear expressions are more complex than other machine learning methods, and the related process is difficult to learn^[48]. Clearly, linear regression is appropriate for selecting the most appropriate descriptor. In this paper, we gained 55%/45% of the optimized data sets, which nicely balanced the accuracies and overfitting of the machine learning model. In the end, the important descriptor was gained by comparing the effectiveness of models for various descriptors.

Performance evaluation

In order to assess the prediction accuracy and model performance, we employed the mean absolute error (MAE), mean square error (MSE) and coefficient of determination (R²). Simply, the smaller MAE and MSE values to 0 and the bigger R² values to 1 suggest the higher prediction accuracy and better model performance. The corresponding equations can be summarized:

Where, n indicates the sample size, indicates experimental value and y_j indicates predicted value, and is the average value.

Extreme feature engineering

Following data pretreatment and feature transformation, the amount and quality of the description dataset must be checked further. Feature processing plays an important role in feature engineering and is also the most time-consuming aspect of data analysis. Because feature processing lacks a defined phase, such as algorithms and models with greater technical knowledge and compromises, there is no unified feature processing way. Fortunately, scikit-learn offers a more comprehensive feature processing approach, which includes data preparation, feature selection, dimension reduction, and so on^[50]. Scikit-learn is a free and open-source machine learning library licensed under the Berkeley Software Distribution license^[51]. Thus, in this paper, the python package scikit-learn was used for data pretreatment, feature transformation, feature processing, machine-learning model training and model performance evaluation^[51-53]. Feature selection is the process of eliminating duplicate and unnecessary characteristics from a data collection, determining the important features in the data set, and eventually obtaining the feature subset^[51]. Wrapper methods are common methods for feature selection^[54,55]. The basic description of wrapper methods is:

• Step 1: A subset of features is chosen to train the model. The model here usually refers to a machine learning algorithm, also called an objective function;

• Step 2: Evaluate the model with a validation dataset;

• Step 3: Perform the above operations on different feature subsets based on some search algorithm;

• Step 4: Based on the evaluation results, the best feature subset is selected.

Clearly, the method for finding the optimal descriptors subset belongs to the family of greedy search algorithms. Wrapper methods include three common selection methods, such as sequential feature selection (SFS)^[56], exhaustive feature selection^[57] and recursive feature elimination^[58]. Among them, SFS includes two algorithms, such as sequential forward feature selection algorithm and sequential backward feature selection algorithm. Sequential forward selection algorithm is about execution of the following steps to search the most appropriate features out of N features to fit in K-features subset. Instead of trying to model all the feature combinations, the sequential forward feature selection algorithm aborts the search by finding an optimal solution, which greatly reduces the computational effort^[56]. Therefore, we adopted the sequential forward feature selection algorithm to perform feature selection. In this work, gradient augmented regression (GBR) was used as the objective function.

The extreme feature engineering created many descriptors, and was followed by a preliminary screening of these descriptors. By analyzing the Pearson correlation coefficient, the top 50 descriptors with the highest Pearson correlation coefficient were successfully selected. Figure 6A shows the Pearson correlation coefficients for different new-constructed descriptors. Figure 6B shows the Pearson correlation coefficients for the as-selected 50 descriptors. Figure 6C shows the Pearson correlation map of the as-selected 50 descriptors and stability. Figure 6D shows the trend between the prediction effect of GBR models and the descriptor number.

Figure 6. (A) Indicators of Pearson correlation coefficients for distinct descriptors with varying sequence numbers and stability. The horizontal axis displays the sequence number for the descriptors while the vertical axis is a reference to the relative Pearson correlation coefficient. (B) Indicators of Pearson correlation coefficients for the selected 50 distinct descriptors with varying sequence numbers and stability. The horizontal axis displays the sequence number for the descriptors while the vertical axis is a reference to the relative Pearson correlation coefficient. (C) Pearson correlation map for the selected 50 descriptors and the stability. The color bar on the right represents the correlation coefficient. (D) R² values of GBR models are used to evaluate machine learning algorithms. There are values of descriptors on the horizontal axis, and R² values for GBR models on the vertical axis.

GBR is an enhancement to the Boosting algorithm^[59]. Boosting is a type of integrated machine-learning algorithm that transforms the poor learner into the strong learner. Each sample was initially allocated an equal weight value in the Boosting algorithm^[60,61]. Because each training produced a significant change in the values of data points, the weight values were processed by adding mis-splitting points at the end of each step, and then N iterations were done to obtain N simple base classifiers^[62,63]. Finally, the N basic classifiers acquired were weighted together to form a final model.

The distinction between GBR and Boosting is that each GBR computation is designed to minimize the last residual. To reduce the associated residuals, a new model must be created in the gradient’s orientation to reduce the residuals. As a result, in GBR, each new modeling aimed to reduce the previous model residuals in the gradient orientation^[64], the associated loss-function negative gradient was employed as the estimated value of the residual in the GBR algorithm, and then the regression tree was fitted^[65,66]. As a weak classifier, GBR typically employs a fixed size regression tree. The capacity to analyze mixed data and create models with complicated functions are two properties of the regression tree that make it more accurate in the promotion process^[64]. The GBR model is as follows:

Here, r_m is the weight, T(x; θ_m) is the regression tree, θ_m is the parameter of the regression tree and m is the number of trees^[67]. The GBR models were built iteratively in this paper. The best descriptors were chosen based on the coefficients of descriptor significance, and this process was repeated for the other descriptors until all descriptors had been explored. Finally, the optimal descriptor subset consists of 23 most important descriptors, as shown in Table 3.

Table 3

Comparison of three evaluation indicators and brand-new descriptors chosen by the GBR model

Method	No.	Descriptors	Evaluation index
			R2	MAE	MSE
GBR	1		0.716913	0.230453	0.110132
	2		0.713565	0.231627	0.111434
	3		0.715694	0.230980	0.110606
	4		0.712077	0.234185	0.112013
	5		0.711042	0.235473	0.112416
	6		0.699977	0.240292	0.116721
	7		0.707636	0.235522	0.113741
	8		0.707648	0.238117	0.113736
	9		0.687836	0.241530	0.121444
	10		0.711809	0.233404	0.112117
	11		0.703372	0.239687	0.115400
	12		0.714631	0.234160	0.111020
	13		0.673522	0.248478	0.127012
	14		0.645613	0.247656	0.137870
	15		0.691656	0.243092	0.119958
	16		0.691806	0.242466	0.119899
	17		0.669111	0.246630	0.128729
	18		0.708067	0.235591	0.113573
	19		0.708067	0.235591	0.113573
	20		0.703448	0.238094	0.115370
	21		0.585948	0.274807	0.161083
	22		0.683018	0.244765	0.123318
	23		0.683018	0.244765	0.123318

GBR: Gradient augmented regression; MAE: mean absolute error; MSE: mean square error.

Figure 6D depicts the relationship between GBR model prediction effect and descriptor numbers. It is clear that when the number of descriptors was raised, the prediction effect of the GBR models grew and eventually stabilized^[68,69]. Clearly, the best effect of the GBR models was obtained when the optimum subset of 23 descriptors was employed. Figure 6B depicts the Pearson correlation coefficients for the 50 descriptors chosen and the stability, which ranges from 0.845766 to 0.839595 with minor variations. These descriptors are strongly correlated with the stability of perovskite compounds. The key to understanding the structure-composition-property relationship is to choose the best relevant description.

Following the selection of 23 key descriptors through SFS, the ideal subset of descriptors was chosen based on the three evaluation indices of the GBR model to train the linear regression model, as shown in Table 3. After a large number of experiments, the results showed that these fluctuations were within the range of 3%. It is apparent that the descriptor exhibited the highest R² of 0.716913, the lowest MAE of 0.230453 and the lowest MSE of 0.230453, respectively, indicating the best model performance. The R², MSE and MAE values for descriptors , , , , , and are exceptionally small. The other 16 descriptors exhibited poor prediction effects and were not considered in the subsequent work. For these 7 different descriptors, the performance of the model was similar, so we chose 7 descriptors, , , , , , and , to train the linear regression model.

Due to the limited capabilities of experimental and theoretical tools, traditional material discovery has always been a process of trial and error. The widely used tolerance factor (t for short) to measure the stability of perovskite was proposed by Goldschmidt in 1926. t has become a popular descriptor of stability and has accelerated stability screening of perovskite over the past century. It is worth noting that Goldschmidt tolerance factor t has been widely used to predict the stability of perovskite structures based only on a universal formula of ABX₃ with matching ionic sizes of A-site, B-site and X-site^[70]. Its expression is:

Here, r_A is the A-site ionic radius, r_B is the B-site ionic radius, and r_X is the X-site ionic radius. This is a semi-empirical formula with an accuracy of only 70% that gives a rough indication of the stability of perovskite materials. The descriptors constructed in this work are not only related to the A-site ion radius, B-site ion radius, but also related to the lattice parameters, which are considered to be key features related to the stability of perovskite materials.

Automated machine learning

By automated machine learning, we discovered the quantitative relationships between various variables based on a collection of data, which resulted in the construction of a mathematical model and the estimation of unknown parameters. The linear regression algorithm, as an effective machine-learning algorithm, accurately depicted the connection of data via the straight line and was better for fitting expressions in this paper^[71]. Table 4 shows the 7 descriptors , , , , , , and the corresponding evaluation indexes selected by the linear regression model. It is easy to see that the last descriptor of achieved the greatest R² value of 0.735529, the lowest MAE value of 0.224526 and the lowest MSE value of 0.102889. As a result, the descriptor of was chosen as the best descriptor for investigating structure-composition-property relationships in perovskite compounds (ABO₃), which was only related with A-site ion radius, B-site ion radius, lattice constant b, and α angle of the crystal structure.

Table 4

The linear regression model’s specified descriptors and assessment indices

Method	No.	Descriptors	Evaluation index
			R2	MAE	MSE
Linear Regression	1		0.735529	0.224526	0.102889
	2		0.732250	0.226233	0.104165
	3		0.732245	0.226250	0.104167
	4		0.728364	0.228253	0.105677
	5		0.721646	0.231526	0.108291
	6		0.718913	0.109353	0.235988
	7		0.718752	0.235184	0.109416

MAE: Mean absolute error; MSE: mean square error.

The straightforward linear equation is intended to represent the structure-composition-property relationships of ABO₃ perovskite compounds after obtaining the optimal descriptor using a linear regression model. The following is the equivalent formula:

Where d_i is the final descriptors, F is the stability, f(d₁, d₂,…d_n) is a linear representation of the structure-composition-property connection. A simple linear expression was produced using linear regression analysis as follows:

where k = 0.1485 and z = -0.0380 are the coefficient values. Following linear regression fit and comparison with DFT calculation value, as shown in Figure 7, the dependability of the automated-machine-learning stability expression was validated. The results showed that the effects of A-site ion radius, B-site ion radius, lattice constant b, and α angle of the crystal structure are more significant than that of other variables. The equation of showed the relationship between structure, composition and property in perovskite oxides. Our technique produces a more accurate expression than the semi-empirical formula. In a nutshell, the novel approach may be utilized to investigate the structure-composition-property relationships of ABO₃ perovskite oxides.

Figure 7. Scatter plots showing correlations between the best descriptor and stability for ABO₃ perovskite compounds. The blue line showed anticipated stability values, whereas the scatter points reflect actual stability values.