Combining DFT with ML to study size specific interactions between metal clusters and adsorbates

TL;DR

This paper presents a machine learning model combined with DFT that accurately predicts interactions between adsorbates and metal clusters using simple, invariant interatomic distance descriptors, significantly reducing computational costs.

Contribution

The work introduces a transferable ML model using only interatomic distances as descriptors to predict cluster-adsorbate interactions with high accuracy, applicable across different elements and molecules.

Findings

ML model achieves AME ~ 0.05 eV in interaction energy predictions

Model accurately reproduces potential energy surfaces for incoming atoms

Descriptors are invariant to rotation, translation, and permutation

Abstract

To date, density functional theory (DFT) is one of the most accurate and yet practical theory to gain insight about materials properties. Although successful, the computational cost is the main hurdle even today. A way out is combining DFT with machine learning (ML) to reduce the computational cost without compromising accuracy. However, the success of this approach hinges on the correctness of the descriptors. In the present work, we demonstrate that, based on {\it only} interatomic distances as descriptors, our ML model predicts interaction energy between an adsorbate and Al cluster with absolute mean error (AME) $\sim$ 0.05 eV (or less) and reproduces the PES experienced by an incoming atom. Our extensive DFT calculations reveal that atoms experiencing identical environment within a cluster have identical interaction energy patterns. Further, we demonstrate that our model is not specific to Al clusters, and could be applied to clusters of different elements as well. Its application to compute PES experienced by various test atoms and molecules in the vicinity of different clusters proves the transferability of the model not just to clusters of different elements but also to various molecules. The descriptors chosen are invariant to rotation, translation, and permutation yet very simple to compute is one of the most crucial points of the present work.

Full text

1 Introduction

Clusters, being the system of few atoms, have very different properties than their corresponding atomic and macroscopic analogue1, 2, 3, 4, 5, 6. Due to the unique arrangement of atoms, the properties of clusters vary substantially with size. This size sensitivity is the key feature of atomic clusters and is reflected in all their properties like, melting7, 8, 9, 10, their growth pattern11, reactivity towards various molecules12, 13, 14, 15, 16, 17, 18, 19, etc. The reactivity of clusters has been investigated by both experimental and theoretical means20, 21, 22, 23, 24, 25, 26, 27. Roach et al. studied the reaction of anionic Al7-73 clusters with H2O molecule23, 24. They demonstrated that Al12**- cluster adsorbs more than one H2O molecule, while Al13**- does not react with H2O. The electronic closed shell structure (40 electrons) was considered to be the reason for non-reactivity of Al13**-. But, Al23**- which is also a closed shell (70 electrons) does react while Al20**- which is not a closed shell (61 electrons) does not react with H2O. Hence this variation in reactivity has geometric rather than electronic origin . By adding an atom in Al12**-, the reactivity completely diminishes due to the absence of adjacent Lewis acid-Lewis base active sites in Al13**- cluster. In smaller size regime where every atom counts, addition or removal of just an atom dramatically changes its properties. Thus, it becomes very difficult to bring out general trends in this size range. It has been also demonstrated by a few groups that not just size, site also affects the reactivity of clusters 28, 29. And hence, to model the interaction of clusters with an incoming adsorbate, all possible adsorption sites for all the clusters must be scanned, which in turn leads to a prohibitively large number of DFT calculations. To overcome this problem, we have employed data driven algorithms of machine learning to predict the site specific interaction energies for various aluminum clusters. Machine learning is preferred over conventional curve fitting because it brings out the underlying complex relations buried in the data set which are useful in prediction for newer and unseen data points as we will demonstrate in this work.

Use of ML methods in combination with DFT is continuously increasing in the field of materials30, 31. It is being applied for accelerated materials discovery 32, 33, 34, to understand the underlying electronic structure35, to obtain chemical information36, 37, to predict the potential energy functions 38, 39, 40, 41, 42, 43, 44, 45 and so on. Success of any ML model hinges on choosing the right set of descriptors. Descriptors should be such that they can bring out accurate and hidden trends from a data set, and yet be as simple as possible. Thus designing features/descriptors that completely describe the system are very crucial. Lot of efforts are being put in developing fingerprints that systematically relate the structural features of samples to their functional properties in quantitative terms46, 47, 48, 49. These set of features then find varied applications like finding similarity between two structures 50, 51, 52, finding the structure-activity relation for various systems 53, 54, 55, screening the chemical space to discover novel materials of desired properties56, 57 or even predict properties for a given material 58, 59, 60, 61. In a study by Hansen et al., they outlined a number of established machine learning techniques and investigated the influence of the molecular representation on ML methods performance. The best methods achieve prediction errors of 3 kcal/mol (0.13 eV) for the atomization energies of a wide variety of molecules62.

Few groups have recently used ML for in silico design of catalysts and proved the validity of their model against the first principles methods63, 64, 65. Wang et al. recently demonstrated the use of artificial neural networks combined with kinetic analysis for rapid screening of bimetallic catalysts. Through a Machine Learning model, they could capture the underlying complex and non-linear interaction between adsorbate and metal, with reported RMSE $\sim$ 0.2 eV64. Another group, Ma et al. adopted the use of ML to capture interactions of adsorbate on multimetallics for catalyst screening of CO2 electroreduction with an RMSE of $\sim$ 0.1 eV63. Other studies integrating ab initio calculations and ML, for transition metal catalysts screening have reported errors (RMSE) as low as $\sim$ 0.12 eV65.

Owing to the excellent catalytic properties, nanoparticles and atomic clusters have always been objects of interest66, 67, 68, 69, 70. Recently, there is an increasing trend in resorting to ML for discovering correlations between geometric structure and catalytic activities67, 70 of metal surfaces as well as nanoparticles66. Recently an ML scheme was proposed to understand catalytic activities based on local atomic configurations and applied to study direct NO decomposition on RhAu alloy nanoparticles68. A local structural similarity kernel known as a smooth overlap of atomic positions (SOAP) was used to find similarities between two geometries based on structural descriptors. Gasper et al. used the gradient-boosting algorithm, for prediction of CO adsorption energies on Pt clusters69. They built predictive models of site-specific adsorbate binding on realistic, low-symmetry nanostructures, with AME $\sim$ 0.1 eV (with respect to DFT). Descriptors used during the training of the ML model in this study comprised of d-band center energy, s and p band center energies, Bader charges, generalized coordination number, etc.

In the present work, we use both DFT and ML techniques, to predict the interaction energy of Al clusters with H atom as an adsorbate. This interaction is studied as a function of increasing size. DFT investigations bring out the one to one correlation between neighbor distance distribution of atoms in a cluster and their corresponding interaction energy. This strong correlation provides the rationale of choosing distances between adsorbate and the surface atoms of cluster as descriptors to train the ML model. And indeed our model based on ‘only’ distances as descriptors could predict the interaction energy with errors as low as 0.05 eV. Further, the transferability of our ML model is demonstrated by its application to different homogeneous (Na10) as well as bimetallic clusters (Al6Ga6). To validate our model the adsorbate is also replaced by other atoms (N), and molecules (N2, O2, and CO).

2 Computational Details

To determine the overall reactivity of a cluster, site specific interaction energy needs to be evaluated. However, owing to the lack of long range order and highly altered short range order, different atoms within a cluster interact differently with an incoming adsorbate. To quantify this site specific variation, we have computed interaction energy of various atoms like H, N and molecules like N2, O2, and CO with Al clusters.

Also, the interaction of H atom with cluster of another element like Na10 and bimetallic cluster like Al6Ga6 was computed. All these resulted into about 35,000 single point calculations. The adsorbates were placed at on-top position of all the surface sites (i.e. surface atoms) for all the clusters with size ranging from 5 to 80. The GS geometries for all the clusters were taken from previously reported work71, 72, 73, 74. As shown in Fig. 1, the adsorbate was kept along the outward radial vector from center of mass of the cluster to surface atom. The distance of adsorbate was varied between 1.30 Å to 3.00 Å from the surface site. All the calculations were carried out within the Kohn-Sham formulation of DFT. Projector Augmented Wave pseudopotential 75, 76 was used, with Perdew–Burke–Ehrzenhof (PBE) 77 approximation for the exchange-correlation and generalized gradient 78 approximation, as implemented in plane wave pseudo potential based code, VASP79, 80, 81. Cubic simulation cell, with the image in each cell separated by at least 15 Å  of vacuum, was used. Energy convergence criteria of 10*-5* eV was used for SCF calculations.

Data collected from DFT calculations was then used to train a ML model. We used the Gradient Boosting Regression (GBR) algorithm as implemented in the scikit-learn python package82. GBR is a regression technique that uses decision tree based classifiers as weak learners. We used the mean squared error function as our loss function (i.e. the objective function to be optimized). The GBR was selected after comparing it against four other regression algorithms viz. Linear Regression, Ridge Regression, LASSO and Stochastic Gradient Descent (SGD). An exhaustive grid search was carried out to find the best parameter values of an estimator. 5-fold cross validation was performed to test accuracy of the model. AME was used as the scoring parameter during cross validation. Multiple checks like plotting the validation curve and learning curves were used to ensure that the model did not overfit the data.

3 Results and Discussion

In clusters, due to the finiteness of the system, every atom does not experience identical environment unlike atoms in bulk. To quantify this variation, nearest neighbor distribution of every atom in all the clusters was studied. In this distribution, distance (di,j) between every pair of atom i and j in a cluster was calculated.

Fig. 2-(a) shows the geometry of Al13 cluster. The corresponding $13\times 13$ distance matrix for Al13 is shown in Fig. 2-(b). To identify the atoms that experienced identical environment (in terms of neighbor distances), distances were arranged in ascending order in the sorted distance matrix as shown in Fig. 2-(c). A careful look at this matrix revealed that there were only two unique rows. Implying that all the twelve surface atoms in Al13 were grouped in two classes as is also seen in Fig. 2-(d). Similarly, for Al12 cluster, eleven surface atoms were grouped into seven classes as shown in Fig. 2-(e-h). It should be noted that upon addition of just one atom in Al12 cluster, seven different classes merge and form only two classes in the case of Al13, i.e. an asymmetric cluster gets transformed into a highly symmetric one. As evident from the Fig. 2-(d) and Fig. 2-(h), the interatomic distances can be used to indicate how (dis)ordered the cluster is. By ‘ordered’ cluster we mean a cluster with many identical atoms in terms of the chemical environment that they experience. For example, Al13, and Al36 are ‘ordered’ clusters because all the surface atoms are grouped into 2 for Al13 as shown in Fig. 2-(d), and 6 for Al36 classes as shown in Fig. 3-(a). Whereas for all the disordered clusters, more than half of surface atoms experience unique environment like 7 in case of Al12 as shown in Fig. 2-(h) and 23 in Al25 as shown in Fig. 3-(b).

This variation in the arrangement of individual atoms is a characteristic of atomic clusters in this size regime. As noted earlier, Roach et al. demonstrated a substantial change in the reactivity of Al12**- compared to Al13**- cluster towards H2O molecule23, 24. A proof of concept for such experimental studies on reactivity lies in our observation based on a significant change in the symmetry of clusters with addition of just one atom as shown in Fig. 2-(d) and Fig. 2-(h). It brings out the fact that variation in behavior of clusters with changing size originates from their geometries. And hence, motivation of our work lies on understanding the interaction of clusters with incoming adsorbates as a function of changing geometries with size (fixed geometry for a size). In our work, cluster geometries were expressed in terms of the nearest neighbor distribution. We demonstrate that the site specific interaction depends upon the nearest neighbor distribution of that specific atom (or site) within a cluster.

All atoms having identical nearest neighbor distribution within a cluster, interact identically with the incoming adsorbate. To elaborate this point further, in Fig. 4 we show the interaction energies for all the atoms within a cluster, for a few representative sizes along with their nearest neighbor distribution (or interatomic distances) in the inset. In the case of Al13 (see inset of Fig. 4-a) as explained earlier all the surface atoms could be grouped in two classes based on their respective interatomic distances, indicating that an incoming adsorbate would experience only ‘two’ different environments. Further, when interaction energy of these surface atoms with an H atom (as adsorbate) was computed, it was observed that atoms belonging to one class interact identically with the adsorbate, resulting into identical interaction energy as shown in Fig. 4-(a). This is evident from the interaction energy of the adsorbate with all the surface atoms when placed at ‘on top’ position. This one to one correlation between identical nearest neighbor distribution and interaction energy is also observed in larger clusters like Al75 and Al55 as shown in Fig. 4-(c) and  4-(d) respectively.

It has been also observed that two sites result into identical interaction energy pattern, if and only if ‘all’ the interatomic distances are identical. For example, in the case of Al9, for two atoms, their first 6 nearest neighbors are at identical distances and the last two distances differ as shown in Fig. 5. However, it has resulted into two distinct interaction energy patterns for the respective atoms as shown in inset of Fig. 5. This one to one correlation between the nearest neighbor distances and interaction energy pattern is observed in all the Al clusters that we have studied with size ranging from 5 to 80. The same trend was observed when H atom was replaced with N atom, and N2, O2, and CO molecules. Also, when the cluster is replaced with that of other elements like Na10, and for bimetallic cluster like Al6Ga6, this correlation holds. This strong correlation between identical adsorption sites and identical interaction energy can be understood from the perspective of Hohenberg-Kohn’s first theorem. HK’s first theorem describes the one to one correspondence between external potential and charge density and hence energy (a functional of charge density). Identical sites are the ones that have same relative distribution of atoms in the cluster i.e. identical nearest neighbor distribution. Which in turn results into identical external potential when a test atom is placed at appropriate position (as described in computational details) and hence, identical interaction energy.

Establishing a structure-property relation to understand reactivity problems better will greatly reduce DFT based computational efforts. But, coming up with a set of descriptors that have characteristics like transferability, universality, potential to capture accurate trends and yet be simple is still an ongoing area of research. And so, it will be interesting to test if this one-to-one correlation between nearest neighbor distribution and interaction energy could be exploited by employing data driven models at a minimal computational cost. Since in this size regime, properties of clusters vary substantially by addition/removal of just an atom, it becomes important to study the interaction of clusters as a function of size. The size and site specific interaction energy data generated through our extensive DFT computations was used to train the ML model. The data was fed to GBR to predict the interaction energies for all unique adsorption sites of the clusters between 5 to 20 and few selected larger clusters (25, 36, 42, 55, 67, and 75). This is a balanced mixture of ordered and disordered clusters. Taking a hint from the DFT investigations, descriptors that captured this structure-property relation were designed. While modelling the interaction of clusters with adsorbate the (dis)similarity between two adsorption sites had to be captured. And hence, nearest neighbor distribution as seen by the adsorbate was the logical choice of descriptors. The chosen set of descriptors did not represent any homometric pairs as the cluster geometries were fixed, while only the distance between adsorbate and cluster varied. For homogeneous clusters, only distances were used as descriptors while nuclear charge was also included for bimetallic cluster.

A trend of reducing prediction errors with increasing system representation was seen for all the clusters that we had studied. For any surface site of an n atom cluster, there will be n distances as descriptors. These descriptors were distances from adsorbate to all the atoms in a cluster and hence invariant to rotation, translation, and permutation of the system. Since, the cluster geometries are always fixed for any given size, n distances are enough to represent the entire system and hence the external potential. The model was trained each time by gradually including more descriptors i.e. distances. In Tab. 1 we list the variation in AME as a function of increasing number of descriptors for smaller cluster sizes. The variation in AME is correlated with interatomic distances from the H atom. We will discuss this further by closely analyzing the specific case of Al13.

In Fig. 6-(a) we have plotted the AME as a function of number of descriptors (number of NN distances from the surface site) used to fit the model for predicting interaction energies for Al13. As noted earlier, Al13 has only two types of atoms. The difference between these two types of atoms in their nearest neighbor distribution is picked up in the ML model. And hence we observed improvement in AME at distances where these two groups differ from each other i.e. AME reduced from 0.13 eV to 0.08 eV with descriptors up to nn4 versus nn5. A similar jump (decrease) in AME was observed when nn8 was also included as shown in Tab. 1. nn8 is the point at which the two classes further separated. In Fig.  6-(b) ML predicted energies for Al13 are plotted against DFT calculated energies. The AME in this specific case is 0.02 eV.

This correlation between reduction in AME and variation in the nearest neighbor distribution was also observed for Al5, Al7, and Al9 clusters (see Tab. 1). The nearest neighbor distribution and variation in AME plots for the above mentioned cases are shown in Fig. 7. Overall, the machine learning model has picked up the underlying correlation between nearest neighbor distances and interaction energy. As has been discussed, by means of distances as descriptors, we are providing information about the external potential and thus catching the essence of Hohenberg-Kohn’s first theorem.

This correlation between nearest neighbor and variation in AME is strikingly evident and easy to capture in smaller clusters. While it is not so clear in the case of larger cluster due to increased complexity of the systems. And this is what reflects into the AME as a function of descriptors as shown in Tab. 2. It was observed that the variation in AME was inconsistent when descriptors up to nn10 versus all distances (nn100%) were used. For example, in the case of Al36 and Al75, reduction in AME was more than 25% for each of them as shown in Tab.  2. It must be noted that Al36 (see Fig. 3-(a)) and Al75 (see Fig. 4-(c)) are highly symmetric clusters. Whereas for asymmetric clusters like Al25, Al55, and Al67 the reduction in errors were less than 5% as evident from Tab. 2. But for another asymmetric cluster, Al42, the reduction is much larger i.e. about 20% which is similar to that of symmetric clusters. Thus generalization of results becomes difficult for larger clusters. Nonetheless, even for larger clusters the one to one correlation between reduction in AME and increasing system representation still holds. The overall AME reported in our work is $\approx$ 0.05 eV.

Since the line of search for all the results discussed above was restricted along the radial vector, to model a real situation wherein an adsorbate can approach the cluster from any direction, all possible directions had to be scanned. The one to one correlation between identical sites and identical interaction would be difficult to quantify for this situation, as now the adsorbate was not placed only at on-top sites. Nonetheless, the same recipe of descriptors was still legitimate as the distances taken were from the adsorbate to the atoms representing the external potential. And so, the same set of descriptors would capture the change in chemical environment as seen by an incoming adsorbate. Thus, a logical proposition is, the potential energy surface (PES) of an adsorbate in the vicinity of a cluster could be explored with our model. To validate this, we scanned PES at 800 different points (in all possible directions) around a cluster and is shown in Fig 8. We computed the interaction energy of H atom for these randomly selected 800 points on a sphere that enclosed the Al13 cluster at its center (see Fig. 8-(b)). The distance of H atom from the closest surface site of cluster varies between 1.60 Å  to 2.69 Å. This result is particularly important because through this we could predict interaction energy at any point on the PES of the cluster-adsorbate system with AME as low as 0.04 eV.

The success of ML model, in this case, is a proof of concept that nearest neighbor distances are the correct choice of descriptors. It was also found that with only 400 points on the sphere we could achieve same level of accuracy as with 800 points as shown in learning curve Fig. 9. Further as seen from the figure, our model has picked up variation in PES quite faithfully. The minima on the PES represents on-top position which is the most favorable position for the H atom on this surface. Whereas the maxima is the least favorable position and turns out to be a bridge position for the H atom.

To further validate our model we tested it on other clusters like Na10 and Al6Ga6. It was observed that for H atom adsorbed on a highly asymmetric Na10 cluster, our ML model with the same recipe of descriptors predicted the IE with AME $\approx$ 0.038 eV. To demonstrate the universality of our work, calculations performed with different adsorbing species on Al clusters are noted below. When a single N atom was adsorbed at on-top positions on Al13 cluster, not only the one to one correlation was observed again but also the AME from ML model was 0.06 eV, i.e. in same range as our previous results. Further the errors for prediction of IE using the same ML model when molecules like N2, O2, and CO were adsorbed on Al12 turned out to be, 0.045 eV, 0.049 eV, and 0.042 eV respectively. Finally, we also tested the validity of our ML model on bimetallic cluster Al6Ga6. AME for an H atom adsorbed on top of this cluster turned out to be 0.09 eV. This error was obtained based on only structural representation of the cluster. With the inclusion of nuclear charge/ionic radii of both the elements of the cluster, viz. Al and Ga, in the descriptor set, the AME got down to 0.058 eV.

In a nutshell, the descriptors used to train the ML model are as simple as distances between adsorbate and the surface atoms correlating the structure and activity between the two. Successful prediction of interaction energy by means of descriptors that systematically represent the external potential catches the essence of Hohenberg-Kohn’s first theorem. Our approach differs from the previous work in a few key ways: 1. Descriptors chosen were such so as to model the interaction of a cluster and an incoming adsorbate. And hence, the chemical environment that an adsorbate experienced was explored, 2. The descriptors used did not represent any homometric pairs as only unique adsorption sites for fixed cluster geometries were used and finally 3. The ML model was trained on purely the structural representation of the cluster.

4 Conclusion

To summarize we have combined DFT with ML to understand the interaction between an adsorbate and clusters. The key results of the present work are as follows:

1. Our extensive DFT calculations establish a one-to-one correlation between the nearest neighbor distances and interaction energy for small Al clusters. The results, as demonstrated, are generic and applicable to all the clusters and different adsorbates.
2. We employ the GBR model to predict the site specific interaction energies by using ‘only’ interatomic distances as descriptors. The absolute mean errors are about 0.05 eV. With this, we also demonstrate that our ML algorithm picks up the one-to-one correlation between the nearest neighbor distribution and the site specific interaction energies and hence essence of Hohenberg-Kohn’s first theorem.
3. We reproduce the PES for a test atom in the vicinity of the cluster by employing our ML model. To get AME about as low as 0.04 eV we require only 400 single point calculations, which demonstrates that we could circumvent the compute intensive DFT by employing this model. 4. Our descriptors are the interatomic distances, and hence the computational cost is negligible. In conclusion, we designed a set of descriptors that were as simple as nearest neighbor distances and yet the ones that could accurately capture the structure-activity relation between the cluster and adsorbate.

5 Acknowledgements

The authors thank Dr. Leelavati Narlikar for fruitful discussions. CSIR-4PI is gratefully acknowledged for the computational facility. KJ acknowledges DST (EMR/2016/000591) for partial financial support. SM acknowledges UGC for research fellowship. SA acknowledges DST-INSPIRE for research fellowship.