Recent progress in the data-driven discovery of novel photovoltaic materials

Data collection

The goal of the data-driven procedure is to find the underlying correlations between the target properties and input features via ML/DL algorithms. To avoid "garbage in garbage out" in ML/DL practitioners, a reasonable dataset matters more than the data-driven algorithms, suggesting that more efforts are demanded to cleanse and filter the original dataset^[36-38].

Descriptor generation

Descriptors, also known as features, and search space are the two crucial issues in ML. The descriptor types of material structures depend on whether the structures are aperiodic or periodic. Aperiodic structures can be stored in the file types of the simplified molecular-input line-entry system (SMILES) or molecular file (MOL), while periodic ones can be dumped into crystallographic information files (CIFs). Aperiodic system structures, including the pure organic structures of photovoltaic absorbers in OSCs, DSSCs, and metal-organic frameworks (MOF), can be depicted using molecular descriptors, while atomic and crystal structural descriptors are generally used for periodic systems, e.g., the crystal perovskite structures in PSCs.

Data preprocessing

The collected data must be preprocessed to check their consistency and noise, especially for the same experimental data with the same testing conditions reported by different researchers. In addition, string variables that contain not just numbers but also other alphabetic characters (typically represented as categorical variables, e.g., lattice structures of hexagonal, tetragonal and cubic) are usually coded into integers by applying coding algorithms such as the one-hot encoder. If 2 variables are highly correlated, one of them should be removed to reduce the redundancy.

Scaling data is one of the prerequisite steps in data preprocessing to transform the input values into the same range, e.g., 0 to 1, which is optional for the tree-based algorithms that are insensitive to variable ranges. Several common scaling methods are accessible in the Python package scikit-learn (sklearn)^[125]. For example, the standardization method is the most widely used and transforms data to the center with a zero mean and unit variance. Min-max scaling transforms the variables to lie between a given minimum and maximum value, often between 0 and 1.

The third important preprocessing step is to randomly divide a whole dataset into three subsets, usually by arranging a training set for building models, a validating set for evaluation while tuning model hyperparameters, and a test set for final evaluation of the model predicting performance. If the whole dataset is sufficiently large, the three subsets should possess the same distributions as the whole dataset^{[126, 127]}. However, the random splitting method does not operate well in relatively smaller and/or sparsely populated data, and hence the trained model tends to misjudge the test samples. Thus, when encountering a small dataset, the K-fold cross validation (K-CV) can be used to replace the validation set. In the K-CV method, the sample sets are randomly divided into K folds, with one of the folds used as the validating set and the rest of the folds acting as the training set. The divided training and validating sets are conducted K-times such that each of the K-fold data is used as a validating set once, and the average performances of the K models are taken as the trained ML model. If the fold number K is equivalent to the number of samples, then each fold contains only one of the samples and the method is deemed LOOCV.

Feature selection

The next crucial step is the feature selection to determine the critical features highly related to the target values and eliminate the redundant variables. Generally, feature selection methods can be classified into three types, namely, filter, wrapper, and embedded methods^[128].

Filter-type methods evaluate variables that only rely on the general characteristics of the dataset and do not involve any ML algorithm, which is advantageous for low computing costs^{[38, 129]}. For instance, the minimum redundancy maximum relevance method (mRMR)^[130-132] selects the optimal features by inspecting the relevance between the features and the target, and the redundancy among features. The maximum relevance is revealed by searching for the features that have the largest mutual information to the target $$ \mathit{\boldsymbol{y}} $$, which satisfies the following equation:

where $$ \mathit{\boldsymbol{S}}_m $$ represents the selected set comprising $$ m $$ features $$ \{\mathit{\boldsymbol{x}}_{i}, i=1, \ldots, m\} $$ and $$ I(\mathit{\boldsymbol{x}}_i, \mathit{\boldsymbol{y}}) $$ evaluates the mutual information between $$ \mathit{\boldsymbol{x}}_i $$ and $$ \mathit{\boldsymbol{y}} $$ as follows:

The redundant information among the features in $$ \mathit{\boldsymbol{S}}_m $$ is restrained by minimizing their mutual information as follows:

where $$ I(\mathit{\boldsymbol{x}}_{i}, \mathit{\boldsymbol{x}}_{j} ) $$ evaluates the mutual information between $$ \mathit{\boldsymbol{x}}_{i} $$ and $$ \mathit{\boldsymbol{x}}_{j} $$:

Therefore, we can combine Equations (1) and (3) and consider the following simplest form to optimize them simultaneously:

Using Equation (5), Gallego et al.^[132] gave one of the simplest algorithms as follows:

(1) Select one feature.

(2) Calculate its mutual information with the target as the relevance.

(3) Calculate its mean mutual information with other features as the redundance.

(4) Determine the difference between the relevance and redundancy as the mRMR score.

(5) Rank the features based on score.

After ranking the features by mRMR scores, one might propose a threshold and select features where the mRMR scores are higher than the threshold. Furthermore, mRMR scores can be combined with an ML model to select features. For example, based on the features ranked by the mRMR score, the recursive feature addition (RFA)^[133] procedure can be used to determine the best feature subset by adding or removing one or more features, as follows:

(1) Select the top feature in the ranked features.

(2) Train and evaluate an ML model.

(3) Select the top two features in the ranked features to evaluate the new model.

(4) Subsequently, select the top three, four, five, and so on features in the ranked features and evaluate the new model, which results in the optimal feature subset with the best model performance.

The opposite recursive feature elimination (RFE) performs the same procedure but starts from the full feature set and eliminates features from the inverse order. By combining the mRMR filter method and the RFA/RFE procedure, an optimal feature subset for the model construction can be obtained, while the mRMR can be alternated by other filter methods, such as the variance threshold, mutual information, and chi-squared test methods^[128].

Differentiating from filter types, wrapper methods select features depending on the model performance of an ML algorithm and typically iteratively repeat two steps: (1) searching for a feature subset; (2) evaluating the model performance with the feature subset, with the best performance corresponding to an optimal feature subset. In the typical representative genetic algorithm (GA)^{[122, 134-137]} method, each feature subset is regarded as a chromosome and has its own fitness that refers to the model performance of a specified algorithm. The superior/inferior chromosomes are retained/discarded, while new chromosomes are regenerated in each iterative step (known as generation) by mutating and crossing over. The detailed procedure of a GA is as follows:

(1) Generate a population composed of chromosomes. Each chromosome represents a feature subset.

(2) Evaluate chromosomes in this population by an ML model with a loss function. The model performance is set as the score for each chromosome.

(3) Deprecate the chromosomes with low scores.

(4) Crossover a randomly selected pair of chromosomes by exchanging their subparts to generate two new chromosomes and supplement the population.

(5) Mutate a randomly selected chromosome (usually $$ < $$ 1%) by slightly altering part of its features to generate one new chromosome and supplement the population.

(6) Repeat steps (2)-(5) until the maximum step is reached.

Another widely used wrapper type is sequential forward selection (SFS) and its backward counterpart (SBS)^[138]. SFS starts with one feature and finds the best feature that can maximize the performance of a model trained by one feature only, where, in contrast to RFA, every one feature is chosen randomly or iteratively. The second feature (every feature is chosen randomly in the rest features) is then added, the model is trained by two features, and the best second feature is selected. The procedure continues until the best performance is found by testing, which finally selects the desired features. SBS follows the same concept but from the full feature set and removes one feature that can maximize the model performance.

Embedded methods perform the feature selection in the process of training ML models and are specific to some ML algorithms that can export feature scores internally, e.g., tree-based algorithms (decision trees, RFs, and so on)^{[129, 139]}. When constructing a decision tree structure, the change in the Gini index caused by each feature is calculated and the features with high influence on the Gini index can be chosen for the selected set to train the decision tree model^[140]. In addition, in the RF algorithm, multiple decision tree models are combined together, and the important features are determined by the average entropies from the sub-trees. In this regard, the features are sorted by Gini entropies in tree models and the feature subset is then selected via the RFA or RFE procedure. For example, Wen et al. employed embedded methods by combining RFE and tree-based models to select the features, which led to up to nine features remaining^[141]. Li et al. adopted the same method to filter the optimal instrumental features for the bandgap of ABO$$ _3 $$ perovskites, in which the model could reach a stable $$ R^2 $$ value of 0.94 in cross validation with the selected 24 features^[101].

The contributions of features to model predictions can be evaluated by SHAP values^{[82, 142]} Given a full feature set $$ T_n $$ composed of $$ n $$ features and a trained ML model (or a fitness function) $$ f $$ that takes a feature set $$ S $$ comprising $$ m (m\leq n) $$ features as inputs and exports a prediction $$ f(S) $$, the Shapley value $$ \psi_i $$ of a feature $$ \mathit{\boldsymbol{x}}_i $$ is hence defined as follows:

where $$ S\in \{ N|\mathit{\boldsymbol{x}}_{i} \} $$ indicates that $$ S $$ will traverse all the feature subset from the total set $$ T _n $$ but exclude the feature $$ \mathit{\boldsymbol{x}}_{i} $$. The sum of the absolute SHAP values over the entire dataset of a feature represents the feature contribution and hence is used to rank the features. In one of our recent works^[143], the SHAP method was used to determine six and four optimal features for the XGBoost and GBM classification models, respectively, rendering an over 85% test accuracy.

ML model construction

We now consider the core stage to select a suitable data-driven model for describing the relationship between the features and properties comprehensively, which also can be regarded as establishing a mapping function with multiple inputs and one or multiple outputs using ML/DL techniques. Benefitting from the decades of efforts taken by scientists in the fields of computer science, mathematics, and other related fields, abundant choices of user-friendly ML/DL algorithms [Table 3] with powerful predictivity have been publicly distributed and widely used, involving the typical tools, such as sklearn^[125], XGBoost^[50], LightGBM^[144], PyTorch^[145], and TensorFlow^[146], which help materials scientists focus on exploiting feature spaces.

ML algorithms

Linear model

Linear regression. One traditional but still widely used algorithm is linear regression (LR), which expects the target value to be a linear combination of the features. A linear model is usually fitted by reducing the residual sum of squares between the observed and approximate target values via the ordinary least square method. In spite of its simplicity, there are still a large number of applications in photovoltaic fields^{[89, 90, 92, 147-149]}. For example, in the works of Kar^{[89, 90, 150-153]}, the LR algorithm was largely employed to fit multiple robust linear models for DSSCs.

Logistic regression classification. Logistic regression classification (LRC) is proposed to complement the classification form of LR by introducing a logistic function to predict the probability of a certain label^[154-156]. Yu et al. built up an LRC model to determine whether a perovskite film exists after post-treatment, which led to a competitive test accuracy of 84% to 86% of the SVM^[157].

Lasso and ridge regression. To reduce the overfitting problem of LR, the L1 regularization penalty is imposed into the calculation for the residual sum to form the lasso regression^{[158, 159]}, while the L2 regularization penalty is expected to obtain the ridge regression (RR) ^[160]. In the work of Li et al., the lasso model, which was fitted to predict the formation energy of hypothetical perovskite materials based on only composition and stoichiometry information, exhibited a 10-fold cross-validation $$ R^2 $$ of 0.75, which was close to the value for RF of 0.80^[161]. Stoddard et al. trained LR, RR, and lasso models to predict the time when the carrier diffusion length of MAPbI$$ _3 $$ dropped to 85% of its initial value with a mean test error of 12.8%^[162]. After applying a kernel trick into the RR algorithm, the kernel ridge regression (KRR) can be formed. Padula et al. employed KRR to perform multiple models based on electronic properties and FPs to predict the device performances, in which the model targeting PCE values had the largest $$ \rho $$ value of 0.68^[163].

Decision trees

The so-called classification and regression tree (CART) algorithms, also known as decision trees, construct a tree-like structure by a binary recursive partitioning procedure capable of processing continuous and categorical features. The data samples are partitioned recursively into the binary nodes in each step (known as depth) by making a decision based on feature attributes until the number reduces to zero or the depth reaches a specified maximum^[164]. Given the naive operating rule, it is simple to understand and interpret CART models by visualizing their tree structures. For example, Paul et al. trained a decision tree model for inorganic-organic hybrid materials to gain more deep insights into the influence of experimental conditions and perovskite properties on the reaction outcomes^[165].

Ensemble methods

Ensemble methods have gained significant popularity in recent years due to their merits of robustness, stability, and generalization^{[38, 166]}, which particularly refers to the tree-based mathematic approaches of assembling multiple CART models to promote the performance over a singular tree model.

Averaging approach. One common case in ensemble methods is the averaging approach, which leverages the outputs from the several parallelly and independently fitted CART models on average. The base models might be trained on different training sets that are sampled from the whole dataset. When the random sample subsets are drawn for each CART model, the algorithm is called pasting^[167]. When the random sample or feature subsets are drawn with replacement, the method is known as bagging^[168] or random subspaces^[169]. If the random sample and feature subsets are both drawn with replacement, the method is entitled RF^[170]. The base models are organized together and make a collective decision, in which the sklearn package provides the basic module "Voting Class" that is convenient to fill any other model rather than the CART model. Takahashi et al. trained a RF model to predict bandgaps of perovskite materials and estimated 9328 candidates, where 11 undiscovered Li-/Na-based structures had an ideal bandgap and formation energy for solar cell applications^[171].

Boosting approach. Another case in ensemble methods is the boosting approach. The critical idea here is to build a series of CART models in sequence, with each model fitted to reduce the whole bias from the former assembled CART models. The final outputs of the boosting model are then determined by the whole sequentially fitted CART models. Adaboost was the first proposed boosting algorithm whose trait is to repeatedly modify the sample weights in each step of building a new CART model, in which the weights of the samples with large predicted errors are enhanced and each new CART model is trained on the reweighted samples^[172]. The most prevailing algorithm under the boosting theory is GBM and its derivatives, also known as gradient boosting trees (GBTs) ^[173]. Rather than modifying sample weights, each CART model in GBM is trained to predict the bias resulted from the whole former models and the finial outputs are the sum of the whole model predictions. The derivatives aim to reduce the computing cost and promote the fitness of GBM. XGBoost is proposed by imposing a second-order Taylor approximation and L2 regularization into the loss function, which can simplify the procedure of building each CART model^[50]. The light gradient boosting machine (LightGBM) introduces gradient-based one-side sampling and exclusive feature bundling to largely reduce the sample numbers and the feature dimensions to lower the computing and memory costs when dealing with gigantic data^[174]. CatBoost adopts an ordering principle to handle the specified cases that contain a large number of categorical features^{[51, 175]}. Sahu et al. employed RF and GBM models to predict the OSC device performance based on 13 material descriptors, giving the similar cross-validation results for PCE values with $$ R^2 $$ values of 0.78 for GBT and 0.76 for RF^[176].

Support vector machine

SVMs, including support vector classification (SVC) and support vector regression (SVR), are also some of the most widely-used algorithms and have become a must-try method because of their robust performance and fast computing efficiency^{[48, 177, 178]}. With a kernel function, SVC finds a separating hyperplane with the maximum margin in a high-dimensional space, while SVR regresses responses and features in a high-dimensional space and tolerates error $$ \varepsilon $$ on each side of the fitting hyperplane. In the work of Wu et al., the SVR model was fitted to predict the unit cell volume of HOIPs for photovoltaic systems, with an $$ R^2 $$ value of 0.989^[179].

Gaussian processes

Gaussian processes (GPs) for ML are developed based on the Bayesian theorem and Gaussian probability distribution^{[180, 181]}. Unlike in other deterministic ML regressions, GP regression utilizes the Gaussian probability distribution to regress data and express regression results in terms of mean and covariance of the maximal posterior distribution. The following are the responses $$ \mathit{\boldsymbol{y}} $$ used in training and the responses $$ \mathit{\boldsymbol{y}}_* $$ to be predicted to form a joint normal distribution of $$ \left(\begin{array}{l} \mathit{\boldsymbol{y}}\\ \mathit{\boldsymbol{y}}_{*} \end{array}\right) $$ with a zero mean and joint covariance:

(7) $$ \begin{equation} p\left(\begin{array}{l} \mathit{\boldsymbol{y}}\\ \mathit{\boldsymbol{y}}_{*} \end{array}\right) \sim N\left(\left(\begin{array}{l} 0 \\ 0 \end{array}\right), \left[\begin{array}{ll} \mathit{\boldsymbol{K}} & \mathit{\boldsymbol{K}} _{*}^{T} \\ \mathit{\boldsymbol{K}} _{*} & \mathit{\boldsymbol{K}} _{* *} \end{array}\right]\right) \end{equation} $$

where $$ \mathit{\boldsymbol{K}} $$ and $$ \mathit{\boldsymbol{K}} _{* *} $$ are the covariance matrixes for data $$ \mathit{\boldsymbol{X}} $$ and $$ \mathit{\boldsymbol{X}} _{*} $$, respectively, and $$ \mathit{\boldsymbol{K}} _{*} $$ is the covariance correlated between data $$ \mathit{\boldsymbol{X}} $$ and $$ \mathit{\boldsymbol{X}} _{*} $$. The joint distribution is expressed by:

(8) $$ \begin{equation} p\left(\begin{array}{c} \mathit{\boldsymbol{y}} \\ \mathit{\boldsymbol{y}} _{*} \end{array}\right)=p\left(\mathit{\boldsymbol{y}}_{*} |\mathit{\boldsymbol{y}} \right) p( \mathit{\boldsymbol{y}} ) \end{equation} $$

Clearly, the predication of $$ p(\mathit{\boldsymbol{y}}_* | \mathit{\boldsymbol{y}}) $$ is also a conditional normal distribution:

(9) $$ \begin{equation} p\left(\mathit{\boldsymbol{y}}_{*} | \mathit{\boldsymbol{y}}\right)=\frac{p\left(\begin{array}{c} \mathit{\boldsymbol{y}} \\ \mathit{\boldsymbol{y}}_{*} \end{array}\right)}{p(\mathit{\boldsymbol{y}})} \end{equation} $$

The prediction follows the normal distribution $$ p(\mathit{\boldsymbol{y}}_*|\mathit{\boldsymbol{y}})\sim N(\mathit{\boldsymbol{K}}_* \mathit{\boldsymbol{K}}^{-1} \mathit{\boldsymbol{y}}, \mathit{\boldsymbol{K}}_{**}-K_* \mathit{\boldsymbol{K}}^{-1} \mathit{\boldsymbol{K}}_*^T ) $$. The covariance matrix $$ \mathit{\boldsymbol{K}} $$ is calculated by a kernel function. For example, the radial basis function in Equation (3) is the common choice to solve non-linear problems:

(10) $$ \begin{equation} k_{i j}\left(\mathit{\boldsymbol{x}}_{i}, \mathit{\boldsymbol{x}}_{j}\right)=\exp \left(-\frac{1}{2}\left|\frac{\mathit{\boldsymbol{x}}_{i}-\mathit{\boldsymbol{x}}_{j}}{l}\right|^{2}\right) \end{equation} $$

where $$ k_{ij} (\mathit{\boldsymbol{x}}_i, \mathit{\boldsymbol{x}}_j ) $$ represents the element value in covariance matrixes and $$ l $$ is a parameter. The GP models were successfully fitted to predict the PCE values for organic photovoltaic materials^{[68, 182]}.

Deep learning

Artificial neural network. The term deep learning (DL) refers to miscellaneous architectures of neural networks. The artificial neural network (ANN) or multi-layer perceptron has the simplest and most understandable structure, whose structural units comprise the fully connected layers (known as dense layers). As shown in Figure 3A, a deep ANN model is composed of multiple dense layers with a conic distribution in layer length. The input data, starting from the input layer, are processed through all the dense layers by multiplying each parametric weight matrix in each layer, which is finally output as the predicted value. A typical example can be found in the work of Li et al., in which the trained ANN model for PSCs achieved the best performance against the other ML models with the highest $$ R^2 $$ value of 0.97 for bandgap predictions and the highest value of 0.80 for PCE predictions^[183].

Figure 3. Architectures of (A) artificial and (B) convolutional networks.

Convolutional neural network. By appending convolution layers in front of the dense structure, as shown in Figure 3B, the convolutional neural network (CNN) is formed to extract spatial features from images, which can be applied in processing characteristic results that are presented as images, such as from XRD, XPS, and so on^[38]. As mentioned in the descriptor section, the GCNN is suitable for convolving spatial structure information from coordinate data and has been applied in predicting the moduli, formation energy, and bandgap of crystal structures from the Materials Project^[52] by Xie et al.^[120].

Recurrent neural network. Other advanced DL architectures may also have significant potential for materials science applications, though there have been few reported publications. For example, when dealing with sequential data, including various spectra data, the recurrent neural network (RNN), long short-term memory (LSTM) network, and gate recurrent unit network can be exerted to train the relevant model^[184].

Generative adversarial network. The generative adversarial network (GAN)^[185-187] is a sophisticated DL architecture consisting of an unsupervised generator model and a supervised discriminator model, which aims to produce promising candidates for inverse design^[184] Specifically, the goal of a generator model is to fit a function $$ p_{model} (\mathit{\boldsymbol{x}}) $$ that approximates the real sample distribution $$ p_{real} (\mathit{\boldsymbol{x}}) $$ with no direct access to real data points. The discriminator model has access to both the real and fake samples (drawn from the generator model), whose purpose is to differentiate from the real or fake. The error via the discriminator model can be used to train both the discriminator and generator models. Given a well-trained GAN model, we may use the generator model to design reliable candidate structures without human intuition. Though the applications of the GAN model are largely restrained because of the difficulties in converging the pair of models and the need for a significant amount of real data samples of high quality^[185-187], we expect that these bottlenecks will be overcome as data samples accumulate and the GAN develops.

Variational autoencoder. The variational autoencoder (VAE)^[188] is a comparable DL architecture to GAN involving an encoder network and a decoder network, whose novelty is to build a so-called latent space to represent the material structures^[63] Crucially, the encoder network maps the material structures (in the SMILE or CIF format) to vectors in a lower-dimensional space known as latent space, which acts to compress the information from the original data into the vector in latent space. Furthermore, the decoder network performs the inverse operations to decompress the vector to its original form. By training both the encoder and decoder networks to process and reproduce the original data, the VAE model is expected to learn the potential features from the real data samples. Benefitting from the continuous and differentiable vectors in latent space, we can extrapolate and construct new reliable material structures by applying direct search engines (e.g., greedy search), since latent space is a continuous vector space.

Evaluation metrics

Before an ML/DL model is trained, as discussed in section "Data preprocessing" (section 2.3), a dataset is usually divided into a training set, validating set, and test set. The prediction of a trained model on the training set is referred to as the training prediction and the relevant metrics are known as the training metrics. The training metrics usually reveal good performance since the samples are already used in training and therefore cannot be an effective indicator for the model performance. Similarly, the validating and test predictions and metrics can be obtained when predicting the samples in the validating and test sets. Good validating metrics are adopted to optimize the hyperparameters, if any. High performance in test metrics signifies excellent predictivity and generalization abilities.

In regression, the determination coefficient ($$ R^2 $$), as seen in the previous examples, is the most critical and common indicator for model performance, in which a higher value signals better model performance. Given an observation $$ y _i $$ of sample $$ i $$, the corresponding prediction $$ \tilde{y}_{i} $$, and the mean value of the observations $$ \bar{y} $$, $$ R^2 $$ is defined as:

(11) $$ \begin{equation} R^{2}=1-\frac{\sum\left(y_{i}-\tilde{y}_{i}\right)^{2}}{\sum\left(y_{i}-\bar{y}\right)^{2}} \end{equation} $$

The Pearson correlation coefficient ($$ \rho $$) is also often employed in regression and in the correlation examination of features. Given the observation (or feature) $$ \mathit{\boldsymbol{y}} $$, prediction (or another feature) $$ \tilde{\mathit{\boldsymbol{y}}} $$, their covariance array $$ cov(\mathit{\boldsymbol{y}}, \tilde{\mathit{\boldsymbol{y}}}) $$ and their standard deviations $$ \sigma $$, $$ \rho $$ is defined as:

(12) $$ \begin{equation} \rho=\frac{{cov}( \mathit{\boldsymbol{y}} , \tilde{ \mathit{\boldsymbol{y}} })}{\sigma_{ \mathit{\boldsymbol{y}} } \sigma_{\tilde{ \mathit{\boldsymbol{y}} }}} \end{equation} $$

The values of $$ \rho $$ range from $$ -1 $$ to $$ 1 $$. The maximum $$ 1 $$ reveals a perfectly positive linear proportion relation and the minimum $$ -1 $$ shows a perfectly negative linear proportion relation.

Various prediction errors can be adopted to signal the model predicting error. For example, the mean absolute error (MAE) is the mean of the absolute difference between each observation and prediction. The mean squared error (MSE) is calculated from the mean of the squared difference. The root mean squaref are the correctly predicted while the off-diagonal elemen error (RMSE) is the root of the MSE. The value ranges of error metrics are largely dependent on the target range.

In the classification task, the total accuracy can be used to indicate the performance of classification models, which is the division between the correctly predicted samples and the whole samples. To gain more detail, the confusion matrix, also known as the error matrix, can be employed, which is an N-square array (N is the categorical number of labels), as shown in Figure 4A. The columns represent the observed labels and the rows indicate the predicted labels (the definitions of the two axes can be swapped). The element in each pixel expresses how many samples belonging to the observed label are estimated as the predicted label. Apparently, the diagonal elements are correctly predicted, while the off-diagonal elements are all incorrect. The accuracies of a specified class can be obtained by dividing the specified diagonal element by the sum of the corresponding row. With regards to the binary classification that contains the positive or negative labels, as shown in Figure 4B, the correctly predicted positive samples are deemed true positive (TP), while the correctly predicted negative, incorrectly predicted positive and incorrectly predicted negative samples are referred as true negative (TN), false positive (FP), and false negative (FN), respectively. Precision is defined as the accuracy of the positive samples, while the recall score refers to the division of TP over the sum of TP and FP, representing the ability of the classification to find the positive samples with the best value of 1. The F1 score is the combination of the precision and recall score, defined as:

(13) $$ \begin{equation} F1=\frac{2 \times (precision \times recall)}{precision + recall} \end{equation} $$

Figure 4. Classification task for (a) multiple and (b) binary labels. TP: True positive; TN : true negative; FP: false positive; FN: false negative.

The F1 score is interpreted as a weighted average of the precision and recall score, whose best value reaches 1 and the worst is 0.

In a clustering task, the Silhouette coefficient (SC) is the most common indicator, which is calculated as^[189]:

(14) $$ \begin{equation} S C=\frac{1}{N} \sum\limits_{i=1}^{N} \frac{b_{i}-a_{i}}{\max \left(a_{i}, b_{i}\right)} \end{equation} $$

where $$ a_i $$ is the mean distance between the sample $$ i $$ and all other points in the same cluster, $$ b_i $$ is similar to parameter $$ a_i $$ but differs in the next nearest cluster and $$ N $$ is the sample number. A higher SC value indicates better performance of a cluster model, namely, the overlapping areas among clusters are close to zero.

ML model applications combined with domain knowledge

Data-driven progress in PSCs

As discussed in the introduction, despite innumerable merits as absorbers in solar cell devices, perovskite materials, especially in the case of HOIPs, still face imperfections regarding scalability, stability, and environmental pollution. Scalability is mostly related to deposition, film formation, and device integration^{[201, 202]}, which are beyond the scope of this review. The remaining two issues are mainly attributable to the unstable structures and the incorporation of Pb in HOIPs, e.g., the mostly used MAPbI$$ _3 $$, formamidinium lead iodide (FAPbI$$ _3 $$), and their derivatives. Most ML/DL studies, accompanying experimental ones, exploit the leading-edge aspects of promoting stability, lowering the fractions of polluting elements as much as possible, or designing new potential material alternatives. The typical ML publications of PSCs are summarized in Table 5.

As light absorbers, the bandgap is one of the most important properties for HOIPs and can act as a simple and initial criterion to rapidly inspect potential candidates. In this context, Saidi et al. established a complex hierarchical convolutional neural network (HCNN) to predict the bandgaps of ABX$$ _3 $$ HOIP structures with the simple inputs of atomic descriptors^[83]. A total of 380 different compositional ABX$$ _3 $$ HOIP structures (expanded to 862 permutations obtained via rearrangements of the tri-halide moiety) were generated by arranging Cs in addition to 18 organic ions at the A site, Pb or Sn at the B site, and three halogens (excluding fluorine) at the X site. Their bandgaps, lattice constants, and octahedral angles were calculated as the relevant concerned ML targets based on DFT, while the structural coordinates extracted from the relaxed structures were treated as the inputs for the ML models. Two convolutional neural network (CNN) models were initially trained to predict the lattice constants and octahedral angles with RMSE values of 0.01 $$ \mathring{\rm{A}} $$ and 40$$ ^{\circ} $$, respectively. Therefore, an assembled HCNN model was formed by piping these two predicted properties as the partial features coupled with the structural coordinates to the third CNN model that was targeted towards bandgaps, which exhibited a low RMSE value of 0.02 eV. Considering the initial inputs and ultimate output, the model of Saidi successfully predicted the bandgap based on only the information of the atomic coordinates rather than any other complicated DFT calculations. Such work may help us to accelerate DFT calculations by predicting DFT properties via ML models. However, it might be more convincible if the constructed HCNN model can be validated by some external samples to show the generalizability of the model available, since the overfitting problem is very common for CNN-like models.

Similar to the work of Saidi, Li et al. explored the chemical space of ABO$$ _3 $$ perovskite materials based on 758 samples distributed over seven kinds of phases with the relevant targeted bandgaps between 0 and 5 eV from the Materials Project as the training dataset and 1641 materials from OQMD as the validation dataset^[101]. The overall workflow of this work can be seen in Figure 5. In total, 66 descriptors were mainly generated from the pymatgen package^[102] in Python (labeled as basic descriptors), plus the so-called bond-valence vector sum (BVVS) descriptors, while four algorithms were used for the ML models, including the SVM, RF, bootstrap aggregating algorithm (Bagging) and GBT. Eight ML models targeting formation energies and bandgaps were established, with $$ R^2 $$ values for the validation dataset of 0.953, 0.949, 0.960, and 0.964 for the formation energies and 0.674, 0.808, 0.790, and 0.822 for the bandgaps, respectively. Similar to the strategy of Saidi, the prediction of formation energies was then treated as the additional input feature, together with the initial basic and BVVS descriptors (totaling 67), to further predict the bandgap, in which the remolded models, via four algorithms, gained $$ R^2 $$ values of 0.734, 0.821, 0.800, and 0.855, respectively. Finally, a feature selection known as "last-place elimination" based on the GBT model was performed, resulting in an optimal set of 26 features and promoted $$ R^2 $$ values of 0.760, 0.813, 0.817, and 0.856, respectively. The contributions of 26 selected features were ranked by the GBT algorithm, and it was found that the electron number of the d orbital played the most important role, followed by the predicted formation energies. The BVVS on the O site also had a significant contribution, which was chosen to characterize the distortion of the BO6 octahedron. A relatively simple GBT model was built to predict the quantum-based bandgaps of ABO$$ _3 $$ perovskite materials, which exhibited the reliable model generalizability on the validation set. The newly introduced descriptor BVVS revealed its distinct promotion for building ML models, which may inspire us to explore more informative descriptors for perovskite materials.

Figure 5. Overall workflow of progressive learning method. The schematic presents the details of a collection of perovskites and the outlines of a progressive learning workflow, including instrumental variable generation, bandgap ($$ E_g $$) prediction, and results analysis. Reproduced with permission from Li et al., J. Mater. Chem. C 8, 3127 (2020). Copyright 2020 Royal Society of Chemistry^[101].

Compared to the discovery of potential materials by predicting suitable bandgaps, Jin et al. established a classification model to directly identify 2D photovoltaic materials^[200]. To perform the classification task, they collected 98 experimentally reported photovoltaic and 98 non-photovoltaic materials, accompanied by the generation of 22 structural descriptors that were evaluated and ranked by their feature importance. The packing factor ($$ P_f $$), average sublattice neighbor count, Mulliken electronegativity minimum value, and average atomic volume played the leading roles in identifying the PV candidates. The backward selection approach was then employed based on the ranked feature importance to exclude the features with minimal influence on photovoltaic properties, leaving 19 features reserved. Several ML algorithms were employed to construct the models, including the GBT, SVM, RF, Adaboost, stochastic gradient descent classifier (SGDC), CART, and LR. GBT gained the best accuracy, recall, and precision scores (all 100%), while the others performed at $$ \sim $$90%. A total of 3011 PV candidates from 187093 unexplored materials in the ICSD were identified by the GBT model. It is noteworthy that the $$ P_f $$ values of these candidates were concentrated between 0.3 and 0.5, and the candidates with the $$ P_f $$ value particularly fixed at 0.33 had a 30% chance of being PV materials. The authors further filtered 26 materials using a criterion known as dimensionality and computed the theoretical PCE-based DFT methods. As a result, three materials, i.e., Sb$$ _2 $$Se$$ _2 $$Te, Sb$$ _2 $$Te$$ _3 $$ and Bi$$ _2 $$Se$$ _3 $$, exhibited the highest theoretical PCEs. Taking the electronic properties of SB$$ _2 $$Se$$ _2 $$Te as an example, it was found that this outstanding performance might be related to the p-p optical transition in Sb$$ _2 $$Se$$ _2 $$Te enabled by the lone-pair s orbitals of Sb and the built-in electric field induced by the asymmetric geometry. Nevertheless, the further application of the constructed GBT model might be constrained due to the small size of the dataset. It, therefore, might be more reliable if the dataset is expanded or the model is validated by unknown samples or stability tests.

Furthermore, Zhao et al. combined a robotic system, ML, and experiments [Figure 6] to assess the photothermal stability of APbI$$ _3 $$ mixed cation perovskites under different aging conditions^[84]. They fabricated over 1400 APbI$$ _3 $$ perovskite samples with 64 compositional combinations by varying the A-site ion (K, Rb, Cs, MA or FA), stoichiometry, coating method (drop or spin coating), aging temperature (60, 85, 100 or 140 ℃), humidity (0% or 10%) and illumination (dark or light). The figure-of-merit, namely, the studied target, was denoted as $$ T_{80} $$, which indicates the time (in hours) required to decay 20% from its initial value. The GBT algorithm was adopted to perform the model construction, with the lower test set RMSE value of 169 compared to the test set RMSEs for LR (651) and RF (527). SHAP combined with the GBT algorithm was used to interpret the feature importance at different aging temperatures. It was found that the over-stoichiometric condition (e.g., K$$ _{0.05} $$FAPbI$$ _{3.05} $$) led to worse stability caused by the higher defect density. The authors also discovered that the incorporation of Cs was beneficial to the stability of perovskites over 100 ℃ but detrimental under 100 ℃, while the doping of MA was overall neutral for stabilizing perovskites and had a positive effect at low temperatures. Subsequently, the authors performed theoretical simulations to compare the energy costs of perovskite decomposition and the activation energies of possible decomposition pathways, resulting in the same conclusion that the incorporation of Cs in FAPbI$$ _3 $$ could increase the crystal formation energy and simultaneously decrease the gas desorption barrier, while MA exerts the opposite effects. Additionally, the authors fabricated MA$$ _\rm{x} $$Cs$$ _{0.15-\rm{x}} $$FA$$ _{0.85} $$PbI$$ _3 $$ with an n-i-p structure of ITO/(SnO$$ _2 $$: PEIE)/(PCBM: PMMA)/MnSO$$ _4 $$/perovskite/PDCBT/Ta-WOx/Au and the device containing the composition MA$$ _{0. 1} $$Cs$$ _{0.05} $$FA$$ _{0.85} $$PbI$$ _3 $$ maintained 90% of the peak PCE value after 1800 h of continuous operation, in which 10 mol.% organic MA and up to 5 mol.% inorganic Cs/Rb might be a promising incorporation strategy to improve the device stability at below 100 ℃. This work excellently applied ML techniques to accelerate the experimental progress by analyzing the influence of vital experimental conditions based on the GBT model and SHAP method.

Figure 6. (A) Crystal structure of a perovskite with multiple cations, including potassium (K$$ ^{+} $$), rubidium (Rb$$ ^{+} $$), cesium (Cs$$ ^{+} $$), methylammonium (MA$$ ^{+} $$) and formamidinium (FA$$ ^{+} $$). (B) Schematic of HTRobot workflow for automatic synthesis and characterization. The red circles in the bottom panel indicate the five different positions tested on each sample. (C) Detailed workflow of high-throughput operation to evaluate perovskite stability. (D) Photograph of HTRobot system, including (1) a robot arm with four pipettes, (2) a camera and humidity meter, (3) a spectrometer to record the absorbance and photoluminescence, (4) 96-well microplates to mix the precursors, (5) a hotplate, (6) a stock solution of PbI$$ _2 $$, FAI, MAI, CsI and so on, (7) a sample stage, (8) pipette tips, (9) a waste container and (10) a heat sealer to optionally fuse microplates with aluminum foil. Ⅰ, Ⅱ, and Ⅲ show a panoramic view of the setup and top views of the film fabrication and solution preparation, respectively. Reproduced with permission from Zhao et al.^[84] Copyright 2021 Springer Nature.

Considering that the addition of an inert capping layer (CL) might be beneficial to the stability of MAPbI$$ _3 $$, Hartono et al. considered 21 organic salts and 2 X-site anions (Cl/Br) to form potential CL candidates in order to identify whether the CL has the ability to enhance the stability of MAPbI$$ _3 $$ and to probe the underlying mechanisms^[81]. For each CL film, 260 samples, along with their 12 processing conditions, were explored under the same aging test conditions (85% relevant humidity, 85 ℃ aging temperature and 0.16 Sun illumination). The authors photographed the samples every 3 min to record the color changes and defined a key descriptor onset for PSC stability as the time intercepts of the rapid color change from black to yellow. The onset was then labeled as the output of the ML models, while the 12 processing conditions were regarded as features coupled with the structural properties generated from the PubChem database, including molecular weight, partition coefficient, indicating the hydrophobicity/hydrophilicity of the molecules, rotatable bond number, complexity, topological polar surface area (TPSA), hydrogen-bond donor number and element numbers for C, H, Br, N and I. In total, six regression algorithms, namely, LR, KNN, RF, GBT, ANN, and SVM, were involved, in which the RF model gained the best RMSE value of 70.8. The RF algorithm was then combined with SHAP to interpret the model result, showing that the number of hydrogen-bond donors and TPSA were the most critical factors in determining the stability.

Motivated by the feature importance ranking, the authors further compared the top-performing CL material, namely, phenyltriethylammonium (PTEA), which had zero values for hydrogen-bond donors and TPSA, with other CLs via the methods of X-ray diffraction (XRD), scanning electron microscopy, grazing-incidence wide-angle-X-ray scattering and Fourier-transform infrared spectroscopy. The XRD data indicated that a new perovskite phase, (PTEA)$$ _2 $$(MA)$$ _3 $$Pb$$ _4 $$I$$ _{13} $$, was formed on the top film of MAPbI$$ _3 $$. The other results revealed that the top-performing CL stabilized the MAPbI$$ _3 $$ perovskite by modifying the surface structure, coinciding with a suppression in the loss of methylammonium and the formation of both PbI$$ _2 $$ and oxygen-containing compounds at the surface of the perovskite. With the feature analysis via the RF model and SHAP method, the authors successfully discovered the vital features and identified PTEA as the most promising CL material. Combined with the results of the characterization, the new perovskite phase was recognized as the main factor influencing the device stability. Such work illustrates that ML technology can help us to find promising materials rapidly and reasonably and even reveal hidden principles.

In addition to materials discovery, ML has also been applied to quantum dynamics to help uncover complex mechanisms, such as charge carrier trapping in perovskites. For example, Zhou et al. employed the KNN algorithm to analyze the calculated results from ab initio nonadiabatic MD and the most important structural factors for the charge carrier dynamics and bandgap of MAPbI$$ _3 $$^[203]. The work started from pristine tetragonal MAPbI$$ _3 $$ with a 48-atom 1 $$ \times $$ 1 $$ \times $$ 1 supercell and a larger 384-atom 2 $$ \times $$ 2 $$ \times $$ 2 supercell. A 9 ps ab initio MD trajectory for both of the two crystal structures was generated with a 1 fs time step. The nonadiabatic coupling (NAC, proportional to the charger carrier relaxation rate) and bandgap were calculated as the targeted properties for the crystal structure in each trajectory. A total of 414 structural and motional descriptors, including bond lengths/motions, bond angles/motions, dihedral angles/motions, crystal lattice motions, relative orientations, and distances, were generated for the 48-atom system, while 5440 descriptors were generated for the 384-atom system. The pairwise mutual information (MI) between each feature and target was estimated based on the KNN algorithm to reflect the correlations. The various angles/motions of I, Pb, and MA (especially the top three angles of I-I-I, I-Pb-I and Pb-I-Pb) shared the majority of the top highest MI values for both targets of the NAC and bandgap, while the crystal lattice motions showed much less importance than its internal bond and angle descriptors. These results reflected three conclusions: 1) the NAC values depended explicitly on nuclear velocity; 2) MA motions had a strong influence on the nonradiative relaxation since the MA motions shared one part of the top highest MI values; 3) the influence on nonradiative relaxation arose from the geometry of the Pb-I sublattice, mainly including I-I-I, I-Pb-I and Pb-I-Pb angles. The work of Zhou tended to explore the key factors that have impacts on the NAC and bandgaps of different MAPbI$$ _3 $$ structures instead of making model predictions directly. Such research might be more exhaustive if more model validations were accomplished.

Data-driven progress in DSSCs

One of the advantages of DSSCs is the mature process of their device fabrication due to decades of their development and optimization in experimental conditions. Most ML/DL efforts so far have focused on accelerating the discovery of new organic dye sensitizers with notable photovoltaic properties that are promising for leading performance in DSSC devices. The typical ML publications of DSSCs are summarized in Table 6.

Most recently, for the purpose of predicting the PCE values for DSSCs, Krishna et al. [Figure 7] prepared the largest (till 2020) dataset composed of over 1200 dyes that could be divided into seven chemical classes to form the corresponding datasets, including 207 phenothiazines, 229 triphenylamines, 35 diphenylamines, 179 carbazoles, 58 coumarins, 281 porphyrins, and 158 indolines, which cover both metal-based and metal-free dye sensitizers^[90]. The dye structures in the seven datasets were depicted by Dragon software version 7^[86] and PaDEL-descriptor software version 2.21^[91] to generate the descriptors based on their 2D dye structures, containing constitutional information, ring counts, connectivity index, functional group counts, atom centered fragments, atom type E-states, 2D atom pairs, molecular properties and extended topochemical atom indices. Each dataset was split into a training set and a test set using either the Kennard-Stone^[204] or modified k-medoid method^[205] in a ratio of 7:3. The descriptor pool was pre-treated to eliminate the intercorrelated descriptors, followed by feature selection using the in-house program "Best Subset selection v2.1 software". Seven descriptor sets were extracted from the feature selection, where 13 descriptors were selected for triphenylamines, 14 for phenothiazines, 13 for indolines, 12 for porphyrins, 5 for coumarins, 11 for carbazoles, and 4 for diphenylamines. For each training set, five statistically acceptable and robust individual models (IMs) were developed.

Figure 7. Schematic representation of the steps involved in the development of QSPR models. Reproduced with permission from Krishna et al.^[90] Copyright 2020 Elsevier.

To enhance the prediction quality of the test set, the authors further used their in-house intelligent consensus predictor tool^[153] to perform "intelligent" selection based on these five multiple PLS models to complement the shortages of any single model in their test set predictions. Therefore, four types of consensus models (CMs) were developed, in which CM0 referred to as the ordinary consensus model, CM1 leveraged the average of the predictions from the qualified IMs, CM2 was the weighted average predictions from the qualified IMs, and CM3 signified the best model from the selected IMs. The CM3 model was the winning model for the triphenylamine, phenothiazine, indoline, and porphyrin datasets with determining coefficients $$ R^2 $$ on the test set of 0.61, 0.73, 0.74, and 0.69, respectively. Furthermore, for the carbazole and diphenylamine datasets, the four CM models had nearly the same performance, with $$ R^2 $$ values of 0.75 and 0.84. However, one model in the IMs shows the highest $$ R^2 $$ value of 0.68 for the case of the coumarin dataset rather than the CMs.

Afterwards, the authors discussed the mechanistic interpretations of all the descriptors obtained from the IMs for each dataset. For example, in total, ten descriptors appeared in the five IMs of the triphenylamine dataset, involving NdsN, B06[C-O], B07[O-S], B09[C-S], B06[O-S], C-038, C-043, nN(CO), EAT_Shape_Y, graph density, F05[N-N] and X4Av. NdsN represents the N atom numbers with double and single bonds (=N-), indicating the tendency of the localized $$ \pi $$-$$ \pi $$* transition due to intramolecular charge transfer transition (ICT) from the triphenylamine donor, which showed a negative impact on the PCEs according to the negative variable coefficient in the IM equations. B06[C-O], B07[O-S], B09[C-S] and B06[O-S] denoted the presence/absence of C-O, O-S, C-S, and O-S atom pairs at the topological distances of 6, 7, 9 and 6 respectively. The presence of B06[C-O] led to the bathochromic shift of the absorption spectra and the enhancement of the molar extinction coefficient of the dye. B07[O-S] influenced the reduction of the absorption range and the shortening of rapid $$ \pi $$-conjunction latency. B09[C-S] was related to the delocalization of the $$ \pi $$ electrons and the blue shift of the ICT band. C-038 represented the Al-C(=X)-Al fragment (Al referred to aliphatic groups and X referred to any electronegative atoms like O, N, S, P, Se and halogens), while C-043 represented the X-CR..X fragment (X refers to any group linked through carbon). Both of them had an impact on preventing the back-transfer of electrons from the conduction band of the semiconductor to the redox couple and thus reducing the charge recombination. The descriptor nN(CO) was the number of imides in the dye structures, which would improve the aggregation property of the dye over the TiO$$ _2 $$ surface and promote the recombination reaction between the redox electrolyte and electrons in the TiO$$ _2 $$ nanolayer. The EAT_Shape_Y dealt with size and branching in the molecular structure, which would enhance the bulk of dyes resulting in sensitized wide-bandgap in nanostructured photoelectrode. The graph density indicated the surface area of the dye, which led to the prolongation of the electron injection into the nanostructured TiO$$ _2 $$ from the dye. X4Av and F05[N-N] were the average valence connectivity and the frequency of two nitrogen atoms at topological distance 5, having positive and negative contributions to the PCE, respectively.

For the other six datasets, the authors completed a full analysis of the relationships among the descriptors, dye structures, and PCEs, as detailed in the original article. Inspired by the comprehensive discussions for the seven chemical classes, the authors designed ten coumarin dyes due to their low PCEs compared to all other studied chemical classes, in which the designed dyes showed a 20.68%-43.51% increase in PCE values (8.93%-10.62%) compared with the maximum reported experimental PCE value of 7.4%. Krishna and co-workers carried out a systematic investigation of the relationships between seven common kinds of organic dyes and device PCEs for DSSCs. More than 1200 dyes were collected and divided into seven datasets for building various PLS models. However, most of the PLS models exhibited relatively low $$ R^2 $$ values of 0.61-0.75, except for the one for diphenylamine dyes. The PLS algorithm is widely used to solve linear problems and might not be suitable for these datasets. More non-linear model algorithms, such as XGBoost and SVM, could be considered to enhance the model predictability. Furthermore, the authors evaluated the designed sensitizers via DFT calculations, while the device performance is also largely subject to other factors, such as the material interfaces. Experimental validations are more encouraged to determine the device performance of the designed organics.

Wen et al. not only established an accurate, robust and interpretable ML model for predicting PCEs based on DFT-calculated descriptors, but also performed a virtual screening and the assessment of synthetic accessibility to identify new efficient and synthetically accessible organic dyes for DSSCs^[141]. A database incorporating 223 reported organic dyes with experimental PCEs over 4% was built, along with the relaxed electronic structures optimized at the M06-2X/6-31G(d) level. The input features were comprised of 21 easily obtained descriptors extracted from the ground-state structures and statistical properties, such as orbital levels, atom counts, and dipole moments, which were further augmented by the expensively calculated vibrational, cationic, anionic, and excited-state properties. To achieve a compromise between the calculation costs and model accuracy, two models (models A and B) were built for the next 2 stepwise large-scale screenings, exerting only the simple and all features separately. Four algorithms, namely, RF, GBT, SVM, and ANN, were picked to perform the models. For model A, the $$ \rho $$ values in CV5 were 0.57, 0.57, 0.63, and 0.65, and the $$ \rho $$ values in the test set were 0.68, 0.64, 0.68, and 0.76 for the four algorithms, respectively, signifying that the ANN model had the best accuracy. The extended descriptors were then incorporated to train model B, eventuating $$ \rho $$ values in CV5 of 0.75, 0.76, 0.74, and 0.74, and $$ \rho $$ values in the test set of 0.70, 0.76, 0.77, and 0.78, respectively.

To enhance the prediction accuracy, the heterogeneous ensemble voting regressor model was built from the GBT, SVM, and ANN based on the voting weight of 4:7:4. The GBT-SVM-ANN model achieves partially higher accuracies with $$ \rho $$ values in CV5 and a test set of 0.70 and 0.79 for model A and 0.76 and 0.78 for model B, respectively. Then, 20 donor groups (D), 12 $$ \pi $$ groups ($$ \pi $$), 6 acceptor groups (A), and 6 auxiliary acceptor groups (Aa) were permutated to form 10080 molecular structures in the configurations of both typical electron donor-$$ \pi $$-bridge-electron acceptor (D-$$ \pi $$-A) and D-Aa-$$ \pi $$-A with the additional electron-withdrawing unit (Aa). Among them, 9886 were left with the converged optimizations in DFT calculations. The 2-stage screening were then performed by adopting the two GBT-SVM-ANN models to predict the structures, emanating 500 molecules left with their predicted PCEs of over 8%. In addition to the PCEs, the authors also considered the synthetic accessibility (SA score) developed in reference^[206] based on molecular complexity and finally shortlisted 8 prominent dyes with SA scores of less than 4.

In an earlier study, Venkatraman et al. explored the absorption shift of a dye sensitizer influenced by the adsorption on TiO$$ _2 $$^[207]. A total of 1961 absorption data of dyes adsorbed on TiO$$ _2 $$ in various solutions were collected from $$ \sim $$ 500 studies. The natures of their spectral shift were determined by red shift (R), blue shift (B), and unchanged (N) with the threshold of less than 10 nm of the maximum absorption difference, in which the unchanged (N) was further grouped into NR (positive difference below 10 nm) and NB (negative difference over -10 nm). The authors considered three classification schemes. The first was the B, N, and R classification problem with the class distribution (2:1:1). The second problem involved B and NR with a 1:1 distribution, while the third was the NB and R classification with a distribution of 2.5:1. The atom-bond sequences and topological indices were generated in 2060 numbers as the input features using the ISIDA Fragmentor 2017^[208] and RDKit^[64], leaving 200 features remaining after correlation filtering. The dataset was randomly split into training (75%) and test (25%) sets. Six classification algorithms were implemented, including linear discriminant analysis (LDA), KNN, SVM, CART, RF, and GBT, in which the RF models gained the best accuracies both for the training and test sets: 71% and 76% for the B: N: R classification, 76% and 80% for the NB: R classification and 76% and 80% for the B: NR classification, respectively. In order to test the performance of the ML models, three dyes (quercetin, 2, 5-dihydroxytetraphthalic acid and carminic acid) in 5 solvents, including dimethylformamide, acetonitrile, toluene, tetrahydrofuran, and methanol, were examined (14 cases in total), with 81%, 81% and 71% predicted by the RF models in the B: NR, NB: R and B: N: R classifications, respectively.

Our work concerning the data-driven discovery of novel DSSCs features the ML-aided design of new sensitizer materials based on BODIPY^[93] and N-annulated perylene (N-P)^[92]. Taking the case of BODIPY as an example, we collected a total of 58 BODIPY sensitizers that could be divided into horizontal and vertical types, with both types consisting of 29 samples. In contrast to the work of Krishna^[90], we generated descriptors as much as possible to depict the structures of the sensitizers using Dragon and JChem software, resulting in 5515 dimensions of the features. A GA was employed to filter the descriptors for the two types of dataset, which were used to construct two LR models (horizontal and vertical models). The performance of the two models targeting PCEs achieved correlation coefficients $$ \rho $$ of 0.926 and 0.898 in LOOCV and 0.895 and 0.928 in test validations. It is noteworthy that the feature interpretations were very useful in designing new structures with quantum-based validations.

In the horizontal model, for example, the most important descriptor, Mor14p (see details in Section S2 and Figures S3 and S4), indicated that more conjugated structures and a larger number of C-S pairs contributed to the PCE values, stemming in the attachment of the groups, such as benzodithiophene, dithienothiophene, thiophene and similar. An additional C$$ \equiv $$C bond and methoxy groups that are near the B atom in the BODIPY core were added under the interpretations of the descriptors nTD and F05[O-B]. According to the mapping fragments of the descriptors, new potential sensitizers were then designed based on the sensitizer structure with the highest PCE in the dataset for each type. The designed structures were further validated using quantum-based evaluations, which revealed that the new candidate possessed the more conjugated structures, larger absorption spectra, faster electron injection efficiencies, and better performance in terms of short-circuit current density ($$ J_{\rm{sc}} $$) and open-circuit voltage ($$ V_{\rm{oc}} $$). The model prediction and quantum-based validation of the designed candidates shared the same results regarding the promising performance. Despite the complete model analysis and the continuous DFT validations, two main deficiencies still exist, namely, a lack of sufficient samples and experimental validations, which may constrain the further applications of the models.

The exploitation of new dye structures might have reached a bottleneck due to the scant absorption ability of singular organic molecules. The introduction of a co-sensitizer to expand the absorption capability is a practice alternative to enhance the performance of such devices. Cooper et al. probed the discovery of new co-sensitizer materials with panchromatic optical absorption for DSSCs using a design-to-device approach integrated with a high-throughput screening and text-mining method^[209]. In total, 9431 dye materials were generated via the text-mining software ChemDataExtractor^[210], including their chemical structures, maximum absorption wavelengths, and molar extinction coefficients. A stepwise screening based on statistics was then processed to shortlist the potential dye structures. In the initial stage, small molecules, organometallic dyes, and the materials that have no absorption in the solar spectra were first removed, leaving 3053 organic dyes. Two key structure-property data indicated the presence of a carboxylic acid group and a sufficiently large molecular dipole moment (over 5 D), which were applied to filter the remaining dyes, resulting in 309 dyes being shortlisted. This information suggested that the dyes contain a high-performance DSSC anchoring group, leading to the effective adsorption onto TiO$$ _2 $$ surfaces to create working electrodes, while the latter information was required for the effective intramolecular charge transfer after photoexcitation.

Afterwards, the authors developed a dye matching algorithm for further screening. Based on the known optical absorption peak wavelengths and extinction coefficients, each potential dye combination for co-sensitization could be ranked using a quality score. The algorithm ensured that the dye combination avoided the optical absorption overlap, exhibited panchromatic absorption, and had an improvement compared to any single dye, yielding 33 remaining dyes. Lately, the highest occupied molecular orbital (HOMO) and lowest-unoccupied molecular orbital (LUMO) energy levels were inspected by DFT. The dye candidate pool was reduced to 29 dyes after consideration of the criteria of LUMO energy level greater than -3.74 eV (TiO$$ _2 $$ conduction band in a vacuum) and HOMO energy level below -4.85 eV (I$$ ^{-} $$/I$$ _3 $$ redox potential in a vacuum), which was essential for forming a standard DSSC device integration. At the final screening stage, 5 dyes, comprising C1, 8c, XS6, 15, and H3, were retained, considering the ease of synthesis and availability for the next stage of experimental validations. The PCE ($$ \eta $$) ratio $$ \eta_{\rm{dye}} $$: $$ \eta_{\rm{N719}} $$ was used to indicate the photovoltaic performance of the five potential co-sensitizers compared to a reference sample N719 dye, in which the co-sensitizer combination XS6 and 15 gained the largest ratio of 0.92. The atomic force microscopy (AFM) and X-ray reflectometry were further employed to characterize the co-sensitizers, indicating that the combination XS6 and 15 possessed the lowest aggregate coverage of 0.3%, the smallest dye-layer thickness of 19 $$ \mathring{\rm{A}} $$, and the highest surface coverage over 70%, which correlated to the best performance in the filtered co-sensitizers.

Data-driven progress in OSCs

In spite of the long history of OSC studies, ML-related ones were scarce until 2018, which might be traced to the complex systems whose active layers mostly comprise binary or even ternary organic systems. Benefitting from the widespread of AI techniques and the stringent requirement for more efficient OSCs materials, ML and DL techniques are blooming to accelerate the process of discovering new potential PV materials for OSC devices. The main challenges issued from the AI work in the OSC field are mainly focused on 1) the representations of complex organic structures, particularly in blend systems, 2) the poor performance of the ML models with the currently maximum $$ R^2 $$ in the test set lower than 0.77^[211], and 3) how to apply models to experiments. To date, the blend active layer system, especially for the binary organic framework of the polymer as an electron donor (D) and the non-fullerene acceptor (NFA) as an electron acceptor (A), has achieved the most promising PCE values of over 18% in OSC devices, better than the single or ternary organic system. Most attention in the OSC community has been focused on this organic system. The typical ML publications of OSCs are summarized in Table 7.

Very recently and impressively, Kranthiraja et al. manually collected 566 polymer-NFA organic photovoltaic (OPV) samples from 253 publications before the end of 2018 to predict PCEs^[99]. The descriptors were composed of the materials properties (MP) and FPs of the polymers (p) and NFA (n), in which MP included the HOMO, LUMO, bandgap, and molecular weight. The RF model was built up and examined by CV5 with the highest $$ \rho $$ value of 0.85, compared to the values of 0.59, 0.79, 0.85, 0.84, and 0.81 for ANN, GBT, SVM, KRR, and KNN, respectively. Based on the robust RF model, descriptor importance was calculated for the polymer-NFA OPV materials. It was found that the sum of the importance of polymer-relating MP(p) only accounted for 6.9% in the whole descriptors, leading to a constant $$ \rho $$ value of 0.85 for the RF model after removing the MP(p) descriptors of the polymer, which was encouraged by the good $$ \rho $$ value of 0.78 for the bandgap RF model and 0.73 for the HOMO RF model that was built solely from FP(p). Given the acceptable RF model, the authors performed a virtual screening process suitable for the representative NFA molecules (abbreviated as ITIC and IT-4F) based on 200932 polymers combined from 382 donor units and 526 acceptor units that were furtherly fragmented from the structures of 566 polymer samples, in which only 1098 ($$ \sim $$0.5% in the virtual space) have been reported in current publications.

To corroborate the model predicting result, the second-ranked polymer, labeled as PBDT(SBO)TzH, in the predicted PCE list of polymer-ITIC was selected for the synthesis, which consisted of benzodithiophene as the donor unit and thiazolothiazole (Tz) as that acceptor unit that were solubilized by sulfur-bridged 2-butyloctyl (BO) chains (SBO). However, the experimental PCE values in polymer-ITIC and -IT-4F were only 4.44% and 3.42% compared to the predicted values of 11.1% and 10.5%, which might be traceable from the poor solubility and rapid aggregation that was presumably not considered in the RF model. To ameliorate the flaws of PBDT(SBO)TzH, 4 variants were designed by replacing the SBO group that was responsible for the aggregation behavior and the varying the chains of the alkylthiophene-flanked Tz group that accounted for the poor solubility. One of the designed structures, marked as PBDTTzEH, showed a relatively similar experimental PCE value (10.10%) to the predicted one (11.17%), though the others still exhibited poor experimental PCE values of 2.15%, 3.97%, and 2.34% compared to the predicted values of 10.28%, 10.72%, and 10.70% due to the two unpromoted imperfections.

Experimental characteristics were measured to identify why PBDTTzEH had excellent performance. The $$ J-V $$ curves indicated the highest $$ J_{\rm{sc}} $$ value of 16.47 mA cm$$ ^{-2} $$ among the 5 polymers and the secondarily large fill factor value of 0.65, while the electrodeless Xe-flash time-resolved microwave conductivity test showed the most efficient mobility both for holes and electrons, which certainly correlated with its superior PCE. In particular, according to the images from AFM, PBDTTzEH exhibited a well-interdigitated morphology, which signaled no aggregation. The work of Kranthiraja not only provided a systematic ML analysis for polymer samples but also performed a detailed experimental validation. Such work may establish a good paradigm of how to use ML to accelerate the discovery of new potential polymer materials for the OSC community.

Another similar example could be seen in the work of Wu et al.^[148] [Figure 8], which explored potential donor and acceptor materials for OSCs. They extracted 565 donor-acceptor pairs from 274 publications as the data samples to predict PCEs, in which each structure in both the donors and acceptors was divided into several fragments that were furtherly encoded by FPs. As a result, there were 31, 14, 27, and 14 in number for the 1-4 fragments for the donors and 30, 18, 6, 22, and 35 for the 1-5 fragments for the acceptors. Therefore, the description of each donor-acceptor sample was expressed by the FP's combinations of fragments. For the ML modeling, various algorithms of LR, LRC, RF, ANN, and GBT were performed based on the training data composed of $$ \sim $$85% of all samples, where the RF and GBT models exhibited more satisfactory performance with $$ \rho $$ values of 0.70 and 0.71 than the values below 0.60 for the other models. A 10-fold cross-validation was appended to evaluate the five models, leading to the lowest MAE value of 0.832 for the RF model and 1.653 for the GBT model compared to the values of 2.1-2.6 for the others, which signaled the robustness of the RF and GBT models.

Figure 8. (A) Schematic of collecting experimental data and converting chemical structures to digitized data. (B) Schematic of machine training, prediction, and method evaluation. Reproduced with permission from Wu et al.^[148], npj Comput. Mater. 6 (2020). Copyright 2020 Springer Nature.

By targeting the three high ($$ > $$ 11%)/moderate/low ($$ < $$ 7%) groups divided by the PCE values, the authors represented the classification models, in which the RF model still had the highest accuracy of 60.23%, especially the accuracy of 65.2% in the case of the high-label group. Given the favorable performance of the RF and GBT models, virtual screening was conducted to predict the PCE values of the automatically generated 32076000 donor-acceptor pairs that were permuted from the fragments of the donors and acceptors in the collected sample set. Following the criteria of easy synthesis and highly predicted PCE values of over 10%, the authors selected six donor-acceptor binary systems composed of two experimentally reported donors abbreviated as PM6 and PBDB-Y and three undiscovered acceptors denoted as Y-ThCN, Y-ThCH3, and Y-PhCl. With the exception of the PM6:Y-ThCH3 system (6.67% in experiments vs. 11.14% and 10.41% predicted by the RF and GBT models, respectively), their experimental PCE values are highly consistent with their predicted values, with the experimental PCE values for PM6:Y-ThCN, PM6:Y-PhCl, PBDB-T: Y-ThCN, PBDB-T: Y-ThCH3 and PBDB-T: Y-PhCl of 13.18%, 15.71%, 11.02%, 11.08% and 11.19% along with the predicted values of 11.56%, 13.30%, 12.73%, 12.49% and 12.55% from the RF models and 10.52%, 13.33%, 11.49%, 11.64% and 11.32% from the GBT models, respectively.

As discussed in the descriptor generation section, most input variables in the collected data are relevant to the processes of synthesis or testing, especially when the samples are sourced from experimental publications, which are far away from the molecular structures. To identify the relationship between material structures and their experimental properties, the descriptors depicting active layer structures are essentially needed, while the descriptor-based studies have been involved in the above cases. In such an instance, it is important to generate useful descriptors representing the structural information as the input variables for model building in the whole ML procedure. In the character of an organic structural scheme, particularly a binary or even ternary arrangement, describing such a complex organic absorber system in a numeric language remains a long-term crucial challenge, but also a promising and effective aspect of promoting the performance of ML models considering the large variety of descriptor choices that have been summarized in the descriptor generation section.

In the work of Zhao et al. ^[212], different kinds of descriptors were investigated for their effects in three ML models with the inclusion of KNN, KRR and SVM, in which the authors categorized them into structural (FPs) and physical (quantum-based) properties, including energy levels, molecule size, absorption, dipole moment, rotatable bonds and the partition coefficient between n-octanol and water, which was labeled as the XLOGP3 descriptor. The dataset included 566 organic donor-acceptor pairs composed of 513 donors and 33 acceptors. It is noteworthy that the authors refined the distance definition of the donor-acceptor pairs for the distance concept in the KNN and the kernel expressions in KRR and SVM, based on the linear distance combination weighted by the physical and structural descriptors of donors, acceptors, and whole systems.

Starting from five physical descriptors, including HOMO-D (where D is donor) to predict the PCE values, LUMO-D, LUMO-A (A is acceptor), reorganization energies of the polymer and acceptor, the LOOCV $$ \rho $$ values of KNN, KRR, and SVM were all below 0.5, while the maximum value could reach as high as 0.72 in the case of KNN after adding two structural descriptors. Although more physical descriptors were considered sequentially, appending new variables, even the terms under more computational cost, did not trigger an improvement in the models. This indicates that the ML models solely based on physical properties might not be useful for practical purposes, since most encoded information in the physical properties was already involved in the structural FPs. However, to balance the model interpretation and predictivity that were represented by the physical properties and FPs, respectively, retaining the amount of interpretable and simply calculated physical properties was still necessary in the ML models. Under this consideration, the optimal descriptor selection, though not suitable for all cases of OSCs, might be a mixture of physical properties and FPs.

Laying aside the traditionally and commonly used organic representations, Meftahi et al.^[213] employed the so-called signature descriptors proposed by Pablo et al.^[214] in 2013 into the ML work of predicting the quantum-based properties (such as energy levels) and the Scharber-model-based results^[215] (such as PCE, $$ J_{\rm{sc}} $$ and $$ V_{\rm{oc}} $$) for 344 small molecule and polymer electron donors and acceptors, in which the dataset, named as Harvard Photovoltaic Dataset (HOPV15), was collected by Lopez et al. under their massive DFT calculations^[69]. Complimented by the Cahn-Ingold-Prelog priority rules and directed acyclic graph-based definitions, the signature descriptors could be generated as thousands of substructures for each organic sample regardless of the structural complexity in a matter of minutes. The so-called Bayesian regularized artificial neural network with the Laplacian prior (BRANNLP) algorithm was used to perform the model building, as well as feature selection based on the embedded ability within L1 regression, which was implemented in the CSIRO-Biomodeller package and TensorFlow^{[216, 217]}. The neural network was composed of only one hidden layer coupled with the input layer containing descriptors and the output layer for predicting the target, contrary to the large and complex framework in DL. All the Scharber model-based results were predicted via signature descriptors, leaving the $$ R^2 $$ values of 0.72 (PCE), 0.65 ($$ V_{\rm{oc}} $$), 0.57 ($$ J_{\rm{sc}} $$), 0.87 (HOMO), 0.94 (LUMO) and 0.83 (bandgap) in the training set and the same sequential values of 0.78, 0.58, 0.60, 0.49, 0.67 and 0.65 in the test set. With the exception of the case of $$ J_{\rm{sc}} $$, all the models exhibited robust performance and admirable predictive ability, which successfully leveraged resource-intensive DFT calculations into the larger regions of materials space at much lower computing costs. Due to the structural complexity, more informative descriptors for organics are urgently required for ML studies, in which Meftahi's work has provided a good guide for us.

Despite the broad research in the OSC field, there is a significant lack of ML-related studies for the ternary system due to its complexity, in which only one publication^[211], to the best of our knowledge, could be searched on the Web of Science. In Lee's work^[211], a dataset of 124 fullerene derivatives-based ternary OSCs samples, regardless of the blend formations such as the composition of either one donor/two acceptors (D: A1:A2) or two donors/one acceptor (D1:D2:A) in the active layer, were constructed from the current literature, along with the theoretical orbital energies of donors, acceptors, and the whole systems. Targeting the PCE, the regression models of RF, GBT, KNN, LR, SVM, and the 99 training samples were undertaken, eventuating in the best performance with LOOCV $$ R^2 $$ of 0.66 and test $$ R^2 $$ of 0.77 for the RF model. The other models exhibited the poor LOOCV $$ R^2 $$ values lower than 0.60 and test $$ R^2 $$ values from 0.25 to 0.73. Additionally, the classifier models of RF, extra tree classifier (ETC), KNN, SVM, ANN and DT algorithms were also performed, rendering the highest LOOCV and test accuracies of 79% and 76%, respectively, for the RF model. Given the outstanding performance, the feature importance of the RF model was extracted, showing the largest contribution of LUMO-D1 (denoted as LUMO of D1) and the second of HOMO-D1. The logical flowchart of one subtree in the RF model was visualized, signifying the key attributes of LUMO-D1, HOMO-D1, and HOMO-A1, which was also coherent with the results shown in other experimental studies according to the authors. Despite lacking further extensions and applications, Lee et al. presented guidance and initiation for the ML researchers that may have an interest in ternary OSC systems^[211].

The analysis of the whole OSC device structure based on a large-scale dataset constituting both single and blend active layers had not been reported before the work of David et al.^[218]. A dataset comprising 1850 device data was prepared, in which most were obtained from the Danish Technical University ranging from 2011 to 2017 and the remaining were manually scraped between 2017 and 2019. Regarding device stability, the numeric data of $$ T_{80} $$ that was defined as the time taken for the device to reach 80% of the initial efficiency ($$ E_0 $$) were extracted along with $$ T_{\rm{S80}} $$ that was the time taken for the device to reach 80% of the stabilized value. Fully 17 categorical features were acquired, covering the device structures and materials, encapsulation, substrate type, test protocols, environmental conditions, light sources, and the measurement conditions, such as temperature, light level, bias condition, and relative humidity.

Considering the 2 different testing conditions based on the International Summit on Organic Photovoltaic Stability (ISOS) protocols^[219], the dataset involving 1149 samples after carefully data cleansing were treated in three modes: the full dataset with 1149 samples, the data (155) conducted with light soaking (ISOS-L) that relates to photostability, and the data (489) conducted with dark storage studies (ISOS-D) that provides information on the tolerance of the solar cells to oxygen, moisture, other aggressive atmospheric components naturally in air. The sequential minimal optimization regression (SMOreg) algorithm^[220] was introduced into ML model building since it had the ability to produce the weights, namely the importance, of each feature, and thence help us to understand the feature significance, in which a positive SMOreg weighting corresponds to a positive influence on stability and the vice versa. The SMOreg model based on whole data signaled the LOOCV $$ \rho $$ of 0.739 and the test $$ \rho $$ of 0.713. The models using ISOS-L and ISOS-D data exhibited the relatively higher LOOCV $$ \rho $$ of 0.819 and 0.767 and the test $$ \rho $$ of 0.734 and 0.659.

Given the quantified significance from the weighting in SMOreg, several important features were identified. For instance, the features that most positively influenced the stability, namely, $$ T_{80} $$, in the whole dataset were the choices of the materials in the first transport layer and active layer. Furthermore, the use of LED lights would benefit $$ T_{80} $$, while all the ISOS testing conditions and the light intensity would damage the stability. For the case of ISOS-L, the most influential features were the device components (such as substrate, transport layers and active layers electrodes), layer materials, and light source, while layer materials, device architecture, and encapsulation method were the most affecting attributes in ISOS-D dataset. In addition to the research focus on the stability, the same methodology was also applied to predicting the initial efficiency $$ E_0 $$ using SMOreg and full data (1347 samples), deriving the LOOCV $$ \rho $$ of 0.739 and the MAE of 0.605%, along with the most significant features involving the choice of active layer and tandem configuration.

Identical to the robotic work of Zhao et al. in the PSC field^[84], Du et al. also utilized a high-throughput robot-based platform, "AMANDA Line One", to realize the superior precise control in experimental conditions at a very large scale so as to form high-quality and continuous sample points, which could also be expanded to the optimization in experimental details for any solution organic semiconductor and interface materials^[221]. For this purpose, the authors fabricated around 100 OSC devices within photovoltaic performance (such as PCE, $$ V_{\rm{oc}} $$, $$ J_{\rm{sc}} $$ and FF) and performed 50-h photostability testing varying in ten-dimensional processing parameters covering D (donor PM6): A (acceptor Y6) ratio, concentration, spin speed, active layer, annealing temperature, active layer annealing time, solvent additives, solvent additives volume, electronic transport layer (ETL) materials, ETL annealing temperature and time, which totally consumed only $$ \sim $$70 h. From the point of statistical analysis of the observed data, several optimum processing parameters were exploited, e.g., a D: A weight ratio of 1:1.2, a low thermal annealing temperature and the others for higher efficiency, as well as high spin speed and active layer annealing temperature below 100 ℃ for longer stability. Furthermore, GP was employed to build the ML models to predict the four photovoltaic parameters, obtaining RMSEs in the test set of 1.175 (PCE), 0.012 ($$ V_{\rm{oc}} $$), 0.055 (FF), and 0.903 ($$ J_{\rm{sc}} $$).

Availability of data and materials

Supporting information available: details of outlier detection algorithms; Illustration for descriptors Mor14p, Mor24m, R2s; Atomic parameters.

Authors' contributions

Wrote the manuscript: Lu T

Supervised this manuscript: Li M, Lu W, Zhang TY

Conflicts of interest

The authors declared that there are no competing financial interest.

Financial support and sponsorship

This study was supported by the National Key Research and Development Program of China (No.2018YFB0704400), the Key Research Project of Zhejiang Laboratory (No.2021PE0AC02), and Shanghai Pujiang Program (21PJD024).

Ethical approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Copyright

© The Author(s) 2022.

Supplementary Materials