Formation energy prediction of neutral single-atom impurities in 2D materials using tree-based machine learning

The data of impurity systems is taken from the IMP2D database [48], which has 14 662 samples available with the formation energy; some of which belong to the adsorbate type, while the others to the interstitial type in 44 hosts. These hosts were previously identified from the Computational 2D Materials Database [48–50] as materials experimentally synthesized in 2D monolayer structures. An adsorbate impurity is defined as an impurity adsorbed onto the surface of the host atoms–for a single atom, it is also known as adatom–whereas an interstitial defect is an impurity that absorbs into the host material as an interstitial, as depicted in figure 1.

Zoom In Zoom Out Reset image size

The formation energy ($E^\mathrm{f}$) in IMP2D was computed via the following equation:

$$\begin{align} E^\mathrm{f} = E_\mathrm{tot}^\mathrm{imp} - E_\mathrm{tot}^\mathrm{pristine}-\sum_{i}n_i\mu_i,\end{align} \tag{ 1 }$$

where $E_\mathrm{tot}^\mathrm{imp}$ and $E_\mathrm{tot}^\mathrm{pristine}$ are the total energies of the impurity and pristine (host) structures, respectively, obtained from DFT. n_i and µ_i are respectively the number and chemical potential of the added atom (dopant). In this database, only one atom is added implying $n_{i = 1} = 1$. The Perdew, Burke, and Ernzerhof (PBE) exchange-correlation functional [51], which is the class of generalized gradient approximation (GGA), is used in their DFT calculations [48]. This functional is well known to have limitations in accurately predicting certain properties, such as giving underestimate band-gap energy, and inaccuracies in describing the localized states in transition metals. The formation energy calculation can be obtained more accurately by applying the GGA + Hubbard U [52–54] or a functional containing nonlocal exchange. However, the PBE functional gives advantages in (i) consistency where data calculated on the same level of theory can make direct comparison of results more valid, and (ii) a good balance between computational efficiency and accuracy [48], especially for a large number of calculations.

While IMP2D contains 14 662 samples with formation energy, some of them are filtered out based on the following three criteria. First, the selected samples hosted by TiCO₂ (144 samples) are filtered out because the difference of the host and impurity energies seems unreasonably high. Second, the selected samples need to have the total energy convergence within 10⁻³ [eV]. This was obtained from the parameter conv2 in the database [48]. Third, the impurity with outlier formation energy (1 sample) is filtered out. Then the chemical features as listed in table A1 are extracted for all hosts and dopants. The samples with missing values are removed. As such, only 5906 samples in 43 hosts are left for running the ML models, as listed in table B1 and supplementary section S1 for more details. The procedure is summarized in figure 2. The distribution of the computed formation energy is shown in figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

After the data selection, the next step is preparing features, as shown in figure 2. We make distinction between two sets of features, which are termed chemical and structural features. A chemical feature is a property that involves chemical elements of host and added atoms such as chemical potential, number of atoms, atomic radius and electronegativity, whereas a structural feature refers to one that involves physical position or relative distance between atoms, which need a JL polynomial for preparation.

The structural features are generated by a vanishing Jacobi radial functions [55] in the spirit of Behler–Parrinello radial functions [56]. Using the added atom as a center, the structural feature can be written as a linear combination of vanishing Jacobi polynomials, $\tilde{P}_n^{(\alpha, \beta)}$, as follows:

$$\begin{align} x^{st}_n = \sum_i^{N_{rc}} \tilde{P}^{\left(\alpha, \beta\right)}_n \left( \cos{\frac{\pi r_i}{r_c}} \right),\end{align} \tag{ 2 }$$

where r_i is the distance between the impurity atom and neighbor atom i within the cutoff radius r_c, as shown in figure 4(a). The vanishing Jacobi polynomials are defined as

$$\begin{align} \tilde{P}^{\left(\alpha, \beta\right)}_n\left(z\right) = P^{\left(\alpha, \beta\right)}_n \left( \cos \frac{\pi r}{r_{c}} \right) - P^{\left(\alpha, \beta\right)}_n\left(-1\right),\end{align} \tag{ 3 }$$

in such a way that the functions vanish smoothly at the distance r_c, as shown in figure 4(b). Here $P_n^{(\alpha, \beta)}$ is the Jacobi polynomial of order n for $z \in [-1,1]$. These features are invariant under translation and rotation of the reference frame and they encapsulate structural information around the impurity atom which depends on symmetry, distances, and the number of neighboring atoms.

Zoom In Zoom Out Reset image size

Here, the cutoff radius of interaction is heuristically chosen with the sphere radius of $r_c = 5$ [Å] and $n_\mathrm{max} = 5$. The hyper parameters α and β of the Jacobi polynomials are set equal to zero, which reduces them to the Legendre polynomials. In essence, we have a set of chemical features $\{x^\mathrm{chem}_i\}$ and a set of structural features $\{N_{rc},\ x^{st}_n\}$ to train the ML models.

The ML algorithms explored in this work are based on decision trees, and they include RF [57–59], and gradient boosting frameworks [60, 61] comprising gradient boosting regression (GBR), histogram-based gradient-boosting regression (HGBR) and light gradient-boosting machine (LightGBM). It should be noted that LightGBM is implemented in [62] and the others in Scikit-Learn in [63].

In RF, an ensemble of decision trees is built where each tree makes a decision of output on its own. Then the ensemble of the decision trees takes an average for the final result. In Scikit-Learn [63], each tree in the ensemble can be built from a sample drawn with replacement (i.e. a bootstrap sample) from the training set for a specific size of sub-samples. The splitting at each node during the construction of a tree can be considered with all input features or a random subset of specific size as max_feature. These two sources of randomness are applied to decrease the variance of the forest estimator. In this work, we employed the mean squared errors (MSE) as the cost function to measure the quality of the tree splitting [63]. Additionally, the hyper parameters of the decision tree, such as the maximum depth of the tree, the minimum number of samples in a leaf node, can be tuned to optimize the model.

In the gradient boosting frameworks, which is a generalization of boosting to an arbitrary differentiable loss function, a decision tree is used as a weak learner. The boosting deploys a sequence of weak learners to correct the output. The errors of the previous weak learners are calculated, and another weak learner is built to correct these errors. The algorithms can add some level of stochasticity by sub-sampling of the training data for a decision tree [63]. However, in GBR and LightGBM, the subsampling is drawn without replacement, whereas in HGBR, it is not implemented [63]. The random subset of features is allowed by all methods, leading to stochastic gradient boosting. On the other hand, LightGBM and HGBR make histograms of samples for each feature. This allows the algorithms to leverage integer-based data structures (histograms) instead of relying on sorted continuous values when building the trees. Furthermore, LightGBM is implemented with the gradient-based one-side sampling (GOSS) and exclusive feature bundling (EFB) algorithms which improve a speed for training a large data set. GOSS selects all samples with larger gradients and randomly selects some samples with smaller gradients for training a weak learner. EFB is used to deal with sparse feature space such as one-hot features. It should be noted that GOSS is not applied in this work since a simple stochastic sub-sampling is faster for a small sample size. In all gradient boosting algorithms here, the squared errors as the lost function is used to measure the quality of the tree splitting [62, 63].

After preprocessing data, we split all samples into two sets, namely the dataset used for model construction and test, and the blind-test set; see figure 5. Accordingly, one host can contain many dopant types, where each impurity structure (sample) has only one dopant, but can have multiple configurations. The blind-test set (371 samples) contains only MoS₂ and W₂Se₄ hosts, which have been a prior separated from all other samples. This blind-test set can indeed be different samples. Here, we simply selected two hosts which have different host space-groups. The dataset containing the rest of the hosts, which is used for constructing an ML model, is further randomly split the samples into two subsets, namely the training data for 80% and test data for 20% disregarding the hots. Alternatively, we can employ a stratified data splitting where the percentage of samples for each host is preserved in the training set and the test set. Notably, the blind-test set contains impurity structures which ML models have never learned in training data.

Zoom In Zoom Out Reset image size

ML models used in this work have many hyper parameters, which can be selected to optimize bias and variance of the models. We used training data to optimize these hyper parameters. To that end, samples in the training data are split into five folds for cross-validation to average the root mean square error (RMSE) of each hyper-parameter set. The best hyper parameters are subsequently used to retrain the model by all training data. Finally, the trained models are used to predict the formation energy of impurities in the test data and those in the blind-test set. Three metrics are used to evaluate the regression performance, namely, the coefficient of determination (R²), the mean absolute error (MAE), and the RMSE; see appendix C. These procedures are summarized in figure 5.

In this work, we used chemical and structural features as described in section 2.2. These features will give information to a decision tree for splitting internal nodes. We compare the prediction results in three cases: (i) using the chemical features only, (ii) using both the chemical and structural features, and (iii) using both the chemical and structural features, but the chemical features do not include hostenergy/atom and host-spacegroup.

We expect that the chemical features, such as a chemical potential of a species or electronegativity, contain insights to the formation energy. For a moment, let us consider a chemical potential µ of a species, which is defined as the rate of change of energy with respect to the change of the particle number of the given species. We can interpret the first two terms of equation (1) as the chemical potential of a species in a mixture, i.e. the rate of change of energy (free energy) with respect to the change in the number of atoms that are added to the system. Therefore, the formation energy in equation (1) is obtained by comparing between the chemical potential of a species in a mixture (the first two terms) and chemical potential µ of a species (the third term). Consequently, the formation energy in equation (1) can indicate the structural stability when forming an impurity structure. Since the chemical potential of a species is computationally less expensive to obtain than the formation energy and available in the database, it can be a promising feature.

In similar regards, the electronegativity is a measure of the tendency of an atom to attract a bonding pair of electrons. The features natoms and Z-dopant/total determine the concentration of impurities. Other chemical features in table A1 are also straightforward to obtain and included to add more degree of freedoms for a decision tree to split into nodes. We remark that the hostenergy/atom feature requires DFT calculations, but it is also computationally less expensive than the formation energy because it is just the host property (in this work, it has only 43 hosts). Moreover, we showed in section 3.2 that even though when this feature is removed, the considered ML models still yield comparably accurate prediction; see table 3 in section 3.2.

In figure 6, the energy change from the first two terms of equation (1) is depicted as a function of the dopant chemical potential. If the host and dopant have weak interaction in the sense that their formation energy is small, the energy change should be close to the dopant chemical potential. In contrast, the stronger interaction will make the energy change more deviated from the dopant chemical potential.

Zoom In Zoom Out Reset image size

For all structures in this database, the host energy and the dopant chemical potential are negative. Considering the formation energy for a given host, when a dopant is added to the host and the interaction between them is weak, e.g. an adsorbate type, the total energy of the composite system is combined and often more negative than the host. If the host remains the same, the energy change (the first two terms of equation (1)) should be close and proportional to the dopant chemical potential. Consequently, if the energy change is more negative than the dopant chemical potential, then the formation energy is negative and conversely. As evident in figure S1 in supplementary section S2 for the adsorbate type when the host is fixed, the magnitude of the formation energy is almost linearly decreased as the magnitude of the dopant chemical potential decreases. We observed that the lower magnitude of the dopant chemical potential has higher chance to yield a negative formation energy.

In contrast, if the interaction between the host and the dopant is stronger, e.g. for an interstitial defect, the host is more affected. In this case, there is a higher chance that the composite system has total energy less negative than the host energy. Consequently, the energy change is positive and thus makes the formation energy positive; see figure 6. Even when the host is fixed for the interstitial type, the formation energy versus the dopant chemical potential has relation that tends to be less than that of the adsorbate type; see figure S2 in supplementary section S2. This makes the interstitial type more complicated to predict than the adsorbate type. However, if we plot the formation energy of all hosts together versus the dopant chemical potential, the trends are less pronounced for both types; see figure S3 in supplementary section S2.

For the electronegativity, if a pair is bonded, the covalent bond requires lower relative electronegativity of the atomic pair, while the ionic bond requires a greater one. The stronger bond is expected to have more negative formation energy. It should be noted that, in this database, greater relative electronegativity of the dopant and its neighborhood has higher probability to make a formation. This is consistent with the previous report [64] that the formation energy of an oxygen vacancy in metal oxide materials was predicted with the relative electronegativity and the impurity concentration (electrons of vacated oxygen per total electrons in the structure), demonstrating that these features play essential roles in predicting the formation energy.

Since the features constitute non-linear relationships with the target, and they depend on the clustering of the hosts (e.g. the dopant chemical potential), the tree-based ML methods can deal with the problem better than a linear regression model. When including the structural features, which are constructed from neighboring atoms in the radious r_c around the dopant, then we removed the feature that gives rise to the clustering of the hosts, i.e. hostenergy/atom and host-spacegroup without a significant loss of the prediction score. Note that the linear regression with the ordinary least square gives worse prediction in the same setting with these feature sets; see supplementary section S3.

The samples are divided into three data sets: training data, test data and blind-test data. The ML models' hyper parameters are tuned by 5-folds cross-validation of training data, whose procedure is detailed in section 2.4. The optimized hyper-parameter sets are given in supplementary section S4 for different ML models and feature sets. The prediction results for the adsorbates, the interstitial defects, and both of them with the chemical features alone, or together with the structural features are reported in tables 1 and 2, respectively. Without the structural features, the ML models give lower MAE and RMSE, and higher R² scores for the adsorbate type than the interstitial type, since the interaction of the latter is more complicated as reflected in figure 6. If the chemical and structural features are applied together, the ML models improve their performance for predicting the formation energy for both impurity types and the union of them. Moreover, the models are also able to predict the formation energy for the blind test, which contains the hosts not recognized by the models from the training data.

It is also interesting to note that if only the structural features are used, the predictions are worse than using only the chemical features, since we considered all atoms in the r_c radius as the same species; hence, they cannot differentiate different host and dopant species. However, when we used them together with the chemical features, they improved the ML models' performance. We also remark that the database contains only relaxed structures of impurities. In practice, the structural features should be constructed from the initial structure, or at least the unrelaxed structures should be used in the predictions. We believe that at worst the performance of prediction is as given by table 1 when only the chemical features are used, and it should be improved even if the non-relaxed structures are used for constructing the structural features.

In addition, we trained and tested the aforementioned ML models by the chemical and structural features together, but hostenergy/atom and host-spacegroup are not included; see results in table 3. In this case, the rest of the features are not specific to the host names explicitly. Comparing the results in tables 2 and 3, these ML models still yield comparable prediction accuracy; see also figure 7 where we compared the MAEs of these feature sets when using the LightGBM model.

Zoom In Zoom Out Reset image size

In figure 8, we provide some examples of the prediction of the formation energy by LightGBM plotted versus the computed formation energy from the database. We have chosen to demonstrate our proposed method by LightGBM for its faster training than RF, GBR, HGBR for comparable RMSE. The results from the other algorithms are presented in supplementary section S5. Furthermore, the feature importance, which describes how much each feature is relevant relative to the others and is used to split the internal nodes in the decision trees, is also reported in supplementary section S6.

Zoom In Zoom Out Reset image size

We should point out that the prediction results with stratification–where the dataset and cross-validation are split in such a way that preserves the percentage of samples for each host–are comparable cross-validation scores to those without stratification. Please see supplementary sections S7.

Hence, it can be concluded that a small set of chemical features from fundamental properties of atoms, together with a set of structural features encoded in the JL polynomials, can be successfully utilized to predict the formation energy of neutral single-atom impurities in 2D materials, as indicated by high R² scores of the test and blind test. In this work, we focused our investigation on the adsorbates and interstitial defects because they are available in IMP2D. However, as substitutional defects in 2D materials are of great interests for most applications but they are not available in this database, we applied these models to them separately and included the results in appendix D.

To compare the performance of the selected ML models, we varied the number N of decision trees as $N \in \{10, 20, 50, 100, 200, 500, 700, 1000\}$. Then, the corresponding RMSE score, training time and predicting time are evaluated by the 5-fold cross-validation of training data. Our computing environment has Ubuntu (verison 20.04.3 LTS) operating system with Intel Core i9-7920X CPU (24 cores in total). RF, HGBR, and LightGBM algorithms are run with multi-threading, and the number of threads is fixed to 12, but GBR is only available to run on a single thread. The results for the adsorbates are shown in figure 9(a) for default hyper parameters and in figure 9(b) for optimized hyper parameters; see list of them in supplementary section S4. For these sets of hyper parameters, the training times fall in between a few tenths to ten seconds, and the prediction times are shorter. Even if the training time for a given set of parameters is in order of a few seconds, the model with the shorter training time requires less time-consuming in hyper-parameter optimizations. The default hyper-parameters yield similar RMSE to the optimized ones; however, the latter will give more robustness against different test-data sets. Also notably from figure 9, LightGBM takes the shortest training time for the same RMSE.

Zoom In Zoom Out Reset image size

In addition, the learning curve, which shows a model's performance as the training size varies, is investigated. The results are shown in figure 10 for the selected chemical and structural feature sets; see also supplementary section S8 for other feature sets. We observed that RF yields lower score compared to the other gradient-boosting frameworks (i.e. GBR, HGBR, LightGBM). The interstitial type shows non-vanishing gradient when training samples are increased. Since interstitial type has smaller sample sizes, this is a possible reason for higher RMSE. Hence, the larger training size is expected to improve the prediction performance of the interstitial defects.

Zoom In Zoom Out Reset image size

It is also evident that the formation energy of the adsorbates is more accurately predicted than that of the interstitial defects. We believe it indicates complexity of the formation energy arising from the environment of the added atoms. In the adsorbate case, the interaction of the dopant with the host atoms is relatively weak, so the energy of the composite system is approximately the sum of the energy from the host and the dopant. Hence, for a given host, the energy change after doping is proportional to a dopant chemical potential. In contrast, the interstitial defects experience many influences from neighboring atoms, adding to a more complex situation. As a result, a set of chemical features may not capture the physical nature of such defects sufficiently accurately. Though the addition of JL terms improves the prediction accuracy due to more degrees of freedom in the modeling, it is likely to account for only a fraction of the interaction. Moreover, it seems that the interstitial type needs more training samples than available data.

It is also interesting to note that both impurity types are added by the same impurity into the same host; the difference is just the position of a dopant. Thus, the features that include the structural information can be used to predict them together.

As the distribution of the computed formation energy in figure 3 has indicated, we selected samples in the database which have converged DFT calculations of energy by a convergence threshold 10⁻³ [eV] as mentioned in section 2.1. In the distribution, the formation energy in range $E^\mathrm{f} \gt -5.0$ and $E^\mathrm{f} \lt 15.0$ [eV] includes 38 samples from the total of 5906 samples. Therefore, this region has fewer samples to train the ML models. If we cut off this region, the regression performance of the testing will be improved, especially in terms of MAE and RMSE. This sometimes makes the blind test yield higher score since there might be no sample in the region.

In this work, we account up to the two-body terms and treat all atoms as single species without specifying the different pairs due to the present of multiple species in each host system. The structural features thus can be improved further by multiplying the vanishing-Jacobi polynomials in equation (3) by an appropriate scaling factor, depending on the pair properties before they are summed to construct the feature $x^{st}_n$. Alternatively, summing up to the three-body terms, where the distance and angle between three bodies are taken in to account, is another approach to improve the structural feature.