ARTICLE

https://doi.org/10.1038/s41467-019-12709-1 OPEN

Statistical learning goes beyond the d-band model

providing the thermochemistry of adsorbates on
transition metals
Rodrigo García-Muelas 1* & Núria López 1*
1234567890():,;

The rational design of heterogeneous catalysts relies on the efficient survey of mechanisms
by density functional theory (DFT). However, massive reaction networks cannot be sampled
effectively as they grow exponentially with the size of reactants. Here we present a statistical
principal component analysis and regression applied to the DFT thermochemical data of 71
C1 –C2 species on 12 close-packed metal surfaces. Adsorption is controlled by covalent
(d-band center) and ionic terms (reduction potential), modulated by conjugation and con-
formational contributions. All formation energies can be reproduced from only three key
intermediates (predictors) calculated with DFT. The results agree with accurate experimental
measurements having error bars comparable to those of DFT. The procedure can be
extended to single-atom and near-surface alloys reducing the number of explicit DFT cal-
culation needed by a factor of 20, thus paving the way for a rapid and accurate survey of
whole reaction networks on multimetallic surfaces.

1 Institute
of Chemical Research of Catalonia (ICIQ), The Barcelona Institute of Science and Technology (BIST), Av. Països Catalans 16, 43007 Tarragona,
Spain. *email: [email protected]; [email protected]

NATURE COMMUNICATIONS | (2019)10:4687 | https://doi.org/10.1038/s41467-019-12709-1 | www.nature.com/naturecommunications 1

ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12709-1

H
eterogeneous catalysis holds the key to solving fundamental demonstrate the robustness of the results42. Thus, large FAIR
sustainability issues by introducing renewable compounds databases43 open alternative paths to rationally design hetero-
as a source of chemicals and energy vectors1–3. The geneous catalysts10, by improving existing thermochemical
rational search for new catalysts benefits from the extensive use of models and generating new ones through statistical learning. For
density functional theory (DFT) and kinetic models derived from instance, the formation energies of large molecules in gas phase
it3–10. This procedure requires sampling the reaction network that can be retrieved from neural networks with an accuracy com-
links reactants, intermediates, and products through transition parable to DFT, 0.04 eV MAE44, thus beyond Benson’s model29.
states. For large molecules, such as those involved in biomass An alternative approach is to predict the formation energy of few
valorization processes, the number of intermediates and transition reactivity descriptors from geometric and electronic features, such
states grow exponentially with the molecular size, rendering the as atomic radius, local electronegativity, ionic potential, and the
computational screening of new catalysts impractical. For instance, coordination of the active site9,18,22,45. Particularly, a recent
the decomposition of a medium size C5 –C6 sugar alcohols SISSO study considers features from the adsorbate, the metal, and
encompasses more than 104 species and 105 transition states in its the adsorption site45, providing excellent predictions with one
reaction network11. The positive aspect is that kinetic barriers are DFT evaluation for each metal. However, the method does not
linked to the formation energies of different intermediates4,12. belong to the explanatory class of machine-learning techniques
Thus, the screening may be reduced to obtain the thermo- and thus results are difficult to interpret. A key step to generalize
dynamics for the adsorbed species and couple them with linear- statistical learning models is to extract physical insights from
scaling relationships and microkinetic models to simulate them. For instance, a feature importance analysis18 rediscovered
operation conditions. The upgrade of biomass-derived molecules the differentiated roles of d- and sp-band contributions31,46.
is often done by metals and alloys1,13, and lately attention has Other studies have highlighted the role of electronic and redox
been drawn to the versatile properties of single-atom alloys descriptors for the thermochemistry of transition metal com-
(SAAs) and near-surface alloys (NSAs)9,14–22. The number of plexes47 and oxide-supported single-metal atoms25. Yet, the
combinations is again unlimited and some have shown an potential of statistical learning in heterogeneous catalysis remains
almost continuum of adsorption strengths21. To this end, new largely unexplored9,26.
thermochemical models based on statistical learning may allow a In the present work, we applied principal component analysis
rapid survey of the energies of adsorbed species for faster and Regression (PCA, PCR) on a set of formation energies
screening8,18,23–27. obtained by DFT. From the descriptors so obtained, we retrieved
The pioneering work by Benson established the basis for the d-band center and the redox ability of the metal as the main
thermochemical scaling relationships of gas-phase molecules controllers of the thermochemistry, along conjugation and con-
already in the 60s28. In this formulation, the formation energy for formational effects. With these descriptors and a minimum set of
a hydrocarbon or oxygenated molecule is obtained as the sum of DFT energy evaluations (around two-thousand), we predicted a
the energies stored on C–C, C–O, C–H, and O–H bonds, con- full thermochemical database of 31,000 species adsorbed on
sidering also the contribution from rings, unsaturations, and pure metals, SAAs and NSAs. The methodology reduces the
radicals28. Despite its simplicity, Benson’s model has an number of explicit DFT evaluations by a factor of 20 keeping its
impressive accuracy for small molecules such as hydrocarbons, accuracy. As the procedure is modular it can be adapted or
alcohols, and ethers, the formation energy of which is predicted extended to other systems in heterogeneous catalysis.
with errors lower than 0.05 eV29.
When molecules adsorb on metal surfaces, the interaction has
covalent, ionic, and dispersion contributions. The most studied Results
term is the covalency appearing from the coupling of the metal Interpretation of thermochemical data by PCA. The first step
sp- and d-states with the adsorbate. The sp part depends on the consists in the generation of a well-converged database of for-
species but it is rather constant along the metals. The second one mation energies on late transition metals: Cu, Ag, Au, Ni, Pd, Pt,
gives the metal-to-metal variability and comes from the d-band Rh, Ir, Ru, Os, Zn, and Cd. This was done using the PBE-D2
center and filling30,31. As a consequence, the adsorption energy of functional following the gold standard in DFT42. The formation
a molecular fragment AHx is a linear function of that of its energies, ECx Hy Oz
, are referred to gas-phase reservoirs of
heteroatom, A, and the slope accounts for the valence of AHx 32. methane, hydrogen, and water, Eqs. (1)–(2). In all cases, the
Further linear dependencies have been identified for heteroatoms lowest energy conformation was employed to ensure that the
belonging to the same group in the periodic table, i.e., P* scales PCA includes the information corresponding to conjugation and
with N*33. In addition, the accuracy of the above models can be conformational changes. The data are arranged in a matrix-E in
improved by using site-specific adsorption rules22,32 and the which the rows span the metals i and the columns correspond to
dependence on the local coordination of the adsorption the adsorbates j. As the data matrix is complete, meaning that it
sites22,34,35. In the particular case of very small nanoparticles, the does not have any missing points, it is suitable for PCA. PCA is a
activity modulation is linearly dependent with the local electro- statistical technique that reduces the dimensionality of E by
static potential36. projecting it along the directions of greatest variability, thus
These adsorption energy models have been extended to mul- reducing its noise and aiding to its interpretability48. The process,
tifunctionalized molecules by combining the heteroatom scalings summarized in Fig. 1 and detailed in the Supplementary Methods,
and the Benson model28,37,38. The combination can be centered proceeds as follows: the average adsorption energy for each
either on individual bond energies37 or on the coordination intermediate, μj , is employed to center the adsorption matrix-E
environment of each heteroatom38. Attempts to generalize sim- and get X. We kept the units of X as eV. This matrix is multiplied
plified thermochemical models to other materials are less fre- at the left by its own transpose to get the covariance matrix C,
quent, but for perovskites and transition metal oxides39–41, which is then diagonalized. The eigenvalues of the diagonal
electronic parameters such as the occupancy of eg orbitals and the matrix D are all positive and are placed in decreasing order,
covalency for the oxygen-transition metal bond were deemed consistently with the eigenvector matrix V. Afterwards, V is
descriptors for their catalytic activity. truncated to kmax principal components and multiplied at the left
Still, the reactivity on metals has provided the largest amount by X to get W and T. These matrices contain the descriptors for
of DFT data and benchmarks on the thermochemistry the metals and adsorbates, t ik and wkj , respectively. The

2 NATURE COMMUNICATIONS | (2019)10:4687 | https://doi.org/10.1038/s41467-019-12709-1 | www.nature.com/naturecommunications

NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12709-1 ARTICLE

PCA L1O a
Pt Au
Get thermochemistry data matrix Remove one metal from the
2.0 Pd
E = (Eij); i runs over metals; thermochemistry matrix Ir

ti2 /eV
j runs over species. Ag
0.0 Rh
Os Ni Cu Cd
–2.0 Zn
Ru
Center the data matrix PCR
X = E–μj
Separate the training set into
–10.0 –5.0 0.0 5.0 10.0
predictors and non-predictors ti1/eV
E′ = (Eij ′) ; E″ = (Eij ′)

Get covariance matrix b 0.50

C = Xt X ij
90
0.40 CH2OCH2O*
Apply PCA on E′
D′, V′, W′, T′ O* 60
Diagonalize covariance matrix
D, V 30
0.30
OCH2CHO*
CH3O* 0
Approximate metal descriptors for
0.20 OH*
the validation (prediction) set

w2j /–
Take kmax principal components Tval = Xval W′
get species descriptors {wkj }
W V (k = 1, 2, ..., kmax) 0.10

Approximate species descriptors

{wkj,val} & j,val 0.00
Get metal descriptors {tik} by linear regression on E″ and T′
H* CCHOH*
T = XW Eij ″ = Σk tik′ wkj,val + j,val C*
CCH2*
–0.10
CHOH*

Predict thermochemical data Validate thermochemical data –0.20

Êij = Σk tik wkj + j Êij,val = Σk tik,val wkj,val + j,val 0.00 0.05 0.10 0.15 0.20
w1j /–
Fig. 1 Principal component analysis (PCA) and regression (PCR). The
Fig. 2 Descriptors from principal component analysis. Descriptors for a
formation energy of species j on metal i is obtained from DFT and Eqs. (1)-
metals fðti1 ; ti2 Þg and b adsorbates fðw1j ; w2j Þg, obtained by PCA. The color
(2). These energies are grouped in thermochemistry matrices, in red, and
scale in b measures the robustness of each species of being a predictor (ιj ).
are approximated following PCA or PCR, in green. PCR was validated by
Those marked in brown are more suitable predictors, as they have large
leaving-one-metal-out of the data matrix (L1O). Variables associated with
projections in at least one wjk and a low SD, Supplementary Eq. 12. Those
metals and species are shown in blue and orange, respectively. In black,
marked in yellow are the least suitable. Relevant species are labeled
variables associated with mathematical procedures. A data flow diagram,
including the sizes of all matrices in this study, is shown in Supplementary
Fig. 3 The first metal descriptor, t i1 spans the larger variability, 17 eV,
whereas the second one, t i2 , accounts only for one-third of
adsorption energies can then be retrieved from Eq. (3). this energy span, Fig. 2a. The metals tend to be ordered left-to-


1 right by their position on each period in the table of elements.
xCH4 þ 2x þ y  z H2 þ zH2 Oþ
! Cx Hy O
z ð1Þ Those more resistant to oxidation, such as Pt and Au, appear on
2
the topmost region.


vasp vasp 1 vasp vasp
All adsorbates appear in the right side of Fig. 2b, meaning that
ECx Hy O
z ¼ ECx Hy O
z  xECH4 þ 2x  y þ z EH2  zEH2 O  Evasp

their first weight, w1j , is always positive. The largest terms belong
2
ð2Þ to highly unsaturated carbonaceous species. Therefore, the first
term, t i1 w1j , can be interpreted as the affinity of metal i to form
^ ij ¼ t i1 w1j þ t i2 w2j þ    þ t ik wkj þ    þ t ik wk j þ μj ð3Þ
E covalent bonds with intermediate j. For late transition metals
max max
(groups 8–11), this characteristic can be mapped to the d-band
center, Fig. 3a, which modulates the adsorption strength on such
The accuracy of Eq. (3) is given by the number of principal surfaces30,31. However, the d-band model cannot describe several
components selected; this is, the number of t ik wkj terms. To assess systems that are relevant for catalysis: it does not apply for
the minimum number of terms in the expansion, kmax , we have adsorbates with an almost completely filled valence shell, such as
evaluated two criteria: the MAE and the variance, see Supple- *OH, adsorbing on alloys with almost-filled d-bands46. Besides,
mentary Table 3. In the first case, the MAE stagnates for two the d-band model cannot describe metals from group 12 (Zn, Cd)
terms to a value that falls within the DFT accuracy. At that point, and above, which can be components in high-entropy alloys
98.1% of the variance is already captured. As a result, only two suited for electrocatalysis21. The second weight for the adsorbates,
principal components, kmax ¼ 2, will be used from now on and w2j , is positive for those that bind through an oxygen atom, such
only two descriptors are needed for metals fðt i1 ; t i2 Þg and as O* and *OCH2 CH2 O*, and negative in species that bind by a
adsorbates fðw1j ; w2j Þg, Fig. 2a, b. These descriptors wrap up the *COH center. Therefore, the second term, t i2 w2j , can be
causes of variability in metal-adsorbate bond energies. interpreted as the ionicity of the metal-adsorbate bond. As such,

NATURE COMMUNICATIONS | (2019)10:4687 | https://doi.org/10.1038/s41467-019-12709-1 | www.nature.com/naturecommunications 3

ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12709-1

a b
10.0 Au

Ag Pt
Au 2.0
5.0 Pd
Cu Ir
Ag

ti1/eV

ti2/eV
0.0 10 0.0
Cd Ni Rh
Zn Ni

ti1/eV
0 Cu
Pt Pd
–5.0 Cd
–10 Ir Rh Ru
Zn –2.0
Ru
–10 0 Os
–10.0 εd – εF/eV

–5.0 –4.0 –3.0 –2.0 –1.0 –1.0 0.0 1.0 2.0

d – F/eV Eh /V

Fig. 3 Interpretation of metal descriptors in physical terms. First and second descriptors for the metals, ti1 and ti2 , plotted against a the d-band center and
b the reduction potential50, respectively. The inset shows that d-band center controls the adsorption energy on late transition metals (groups 8–11)30,31,
but it cannot describe the behavior of Zn and Cd. Additional data and plots are provided in Supplementary Table 2 and Supplementary Figs. 1 and 5

t i2 can be mapped to the reduction potential, Fig. 3b. Cross- thermochemistry of a given metal can be estimated from two
correlations do not appear, Supplementary Fig. 1, meaning that principal components obtained from the formation energies of
t i1 w1j and t i2 w2j summarize rather independently the covalent three27 predictors (O*, *OH, and *CCHOH) that capture most of
and ionic contributions, respectively. the variability of the original data matrix. The principal
The fact that 98.1% of the variance is captured by the components define two descriptors for both metals (t i1 , t i2 ) and
aforementioned descriptors, implies that other contributions to molecules (w1j , w2j ).
the thermochemistry, such as conjugation, conformational
changes (different adsorption sites), and dispersion (van der Fast and accurate prediction of thermochemistry via PCR. To
Waals) are already included, as they are related to the two major assess the accuracy of the statistical learning tools, we performed
covalent and ionic terms. For instance, CHCH2 can adsorb as two set of tests, using PCA and PCR-L1O, respectively (Fig. 4). In
monodentate (*CH=CH2 ), tridentate (**CH–*CH2 ), or inter- both cases, we compared the results of the estimated formation
mediate structures. The most stable conformation depends on the ^ Fig. 1) with those obtained by DFT (E,
energies from Eq. (3) (E,
metal affinity to carbon, Supplementary Fig. 4. This means that
Fig. 1). For PCA, the training and prediction sets are equal, as
the adsorbate valence is not necessarily an integer, as it is
they contain the full DFT set: 71 molecules and 12 metals. PCA
normally assumed in heteroatom scaling relationships32. Also,
estimates all the energies within ±0.50 eV and 98% of them lie
small molecules such as *OH can adsorb on fcc, hcp, bridge, and
within ±0.30 eV, with a MAE of 0.08 eV (Fig. 4a, b). The test on
top-tilted sites49, and the preferred site can differ even for
the PCR is stricter as the training set is reduced to only three
chemically similar metals, such as Rh (fcc sites) and Ir (bridge), or
molecules as predictors. We followed a leave-one-out (L1O)
Ru (hpc) and Os (bridge). In other words, the most stable
validation in which we took a subset of 11 metals (training set) to
conformation of the adsorbate is defined by the metal, thus
highlighting the interplay of the metal-adsorbate system. predict the thermochemistry of the 12th one (validation set)
To provide a rapid survey on the adsorption energies of surface (Fig. 1). This matrix of energies is split into two submatrices
species, it is desirable to calculate by DFT only a small subset of containing just the three predictors (E0 ) and the remaining species
intermediates, called predictors. Their number should be at least (E00 ). In total, 792 DFT evaluations are required to get E0 , cor-
equal to the number of principal components kmax 27,48. Choosing responding to 11 metals ´ 71 species plus the empty surface. For
the predictors is not evident from Fig. 2b alone. For instance, the the validation set (Eval ), only four DFT evaluations are needed
simplest set would contain only the heteroatoms C* and O*27, but and correspond to the clean surface and the three predictors. The
EO and EC are mildly codependent. This codependence appears in PCR-L1O starts by applying PCA on the E0 submatrix to obtain
many pairs of adsorbates, as can be seen in Supplementary Fig. 2, T0 and W0 . Then the descriptors for the metals in the validation
although its origin was unknown33,51. Therefore, the predictor set set, t 1;val and t 2;val , are estimated from the DFT formation energies
can be expanded to ensure that the full DFT database, (E in of the three predictors. The descriptors for the remaining species
^
Fig. 1), is properly estimated (E). (w1j;val , w2j;val , and μj;val ) are then found via linear regression of
To select the proper predictor set, we calculated the error Eq. (4) on the training set. Finally, the thermochemistry of the
matrix ðE ^  EÞ that compares the potential energies obtained by validation set is predicted from Eq. (5). This procedure is
DFT to the ones estimated by Eq. (3). Then, the robustness of sequentially run to consider every metal independently as a
each intermediate as a predictor, ιj , was measured following the validation set. Fig. 4c, d compares the results from PCR-L1O with
DFT data, showing that the MAE increases to 0.12 eV and the
procedure detailed in Supplementary Methods. This ιj is indicated
population in the central bars is about 25% smaller than for PCA.
by the color scale in Fig. 2b, in which dark brown indicates at Still, 98% of the estimated energies lie within ±0.40 eV. Thus, the
least one large wkj and a low SD. Therefore, such species are better PCR-L1O methodology only increases the error span by 0.10 eV.
at predicting the thermochemistry of others. Interestingly, we If C* and O* were used as predictors instead of O*, *OH, and
noticed that when the prediction error of O* was positive, the one *CCHOH, the MAE would have rise to 0.16 eV and the max-
of *OH tends to be negative and vice versa. We also chose imum error to ±1.00 eV. Writing the formation energies as linear
*CCHOH as predictor, as it has a pair (w1j ; w2j ) that is almost regression of the d-band center and the oxidation potential
orthogonal to O* and *OH, Fig. 2b. In summary, the full (excluding Os, Cd, and Zn) would increase the MAE to 0.18 eV

4 NATURE COMMUNICATIONS | (2019)10:4687 | https://doi.org/10.1038/s41467-019-12709-1 | www.nature.com/naturecommunications

NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12709-1 ARTICLE

a PCA c PCR-L1O a b
Cu 0.0
6.0 Ag 6.0
Au
–1.0
Ni
4.0 4.0

Ê/eV
Pd

Ê/eV
Eval
Ê/eV

1×3
Pt –2.0
Rh 2.0
2.0 Ir
E Ê E′&E″ Êval Ru 0.0
Os
0.0 11×3 1×68 –5.0
12×71 12×71 Cd
11×68
Zn –5.0 –2.0 –1.0 0.0 0.0 2.0 4.0 6.0
0.0 2.0 4.0 6.0 0.0 2.0 4.0 6.0 EEXP/eV EDFT,val/eV
EDFT/eV EDFT/eV
c d
b d
MAE = 0.08 eV MAE = 0.12 eV 6.0 6.0
300

4.0 4.0

Ê/eV

Ê/eV
200
Count

2.0 2.0

100 0.0 0.0

0 0.0 2.0 4.0 6.0 0.0 2.0 4.0 6.0

–0.5 0.0 0.5 –0.5 0.0 0.5 EDFT,val/eV EDFT,val/eV
PCA/eV PCR-L1O/eV
Fig. 5 Prediction of adsorption energies on metal alloys through PCR. a PCR-
Fig. 4 Accuracy assessment of PCA and PCR. a Error distribution and b L1O estimates of formation energies plotted against experimental data of
cumulative errors for the PCA, taking the DFT data of the 71 molecules on Supplementary Table 4 (orange circles)53. Open circles stand for the
the 12 metals. c, d The corresponding values for PCR-L1O, using only O*, predicted values when coverage effects are not considered. DFT data from
*OH, and *CCHOH as predictors, and leaving-one-metal-out of the pool ref. 52 are provided as benchmark for PBE and BEEF-vdW density functionals,
(L1O). The inset shows the size of the corresponding input (red) and output in dark and light gray, respectively. b Single-atom (SAA), c sub-surface (NSA-
(green) energy matrices SS), and d overlayer (NSA-OL) alloys. Taking the DFT thermochemistry of 71
adsorbates on 12 pure metals and 3 predictors for each alloy (*O, *OH, and
*CCHOH), the adsorption energies for the remaining 68 adsorbates was
and the maximum error to ±0.80 eV. Thus, higher accuracies can
predicted. Then, 1% of these species were randomly selected as validation
be obtained from predictors calculated by DFT than by using
set, calculated by DFT, and benchmarked on b–d
tabulated data.
E00ij ¼ t 0i1 w1j;val þ t 0i2 w2j;val þ    þ t 0ikmax wkmax j;val þ μj;val ð4Þ total of 12
ð15  1Þ ¼ 168 SAA and 336 NSA. For each alloy, the
prediction of the thermochemistry required four DFT evaluations,
^ ij;val ¼ t i1;val w1j;val þ t i2;val w2j;val þ    þ t ik ;val wk j;val þ μj;val corresponding to the clean surface and three predictors: O*, *OH,
E max max and *CCHOH. The alloys whose structures did not converge are
ð5Þ listed in the Supplementary Methods and were removed from the
pool, leaving 165 SAAs and 278 NSAs. These results, collected on
matrix Eval , required 1772 DFT evaluations. With this data, and
The adsorption energies of metals that are better predicted taking the pure metals as a training set (12 ´ 71 = 852 DFT
correspond to Rh and Ir with 0.05 and 0.08 eV MAE, whereas the evaluations distributed between matrices E and E’), a PCR was
worst are Zn and Ag with 0.18 and 0.17 eV MAE. Zn (Ag) has the applied to predict the thermochemistry for the remaining
more exothermic (endothermic) average of formation energies. 68 species on each alloy. As a result, a database of 31,453
Therefore, the binding energy of their adsorbates is normally formation energies was generated (Supplementary File: alloys-
underestimated (overestimated). In general, the predicted values prediction.csv). To benchmark this database, 1% of these
are least accurate for species that have highly endothermic species were randomly selected, calculated by DFT, and compared
formation energies, such as *CHOHCHOH and *CCH. However, with the predictions of PCR. The benchmark, shown in Fig. 5b–d
they would play a minor role when large reaction networks are excluded those species that belong to predictor set. In all cases, the
taken as a whole particularly when different competitive paths are estimates are within the DFT accuracy, with a MAE around 0.19
wrapped through microkinetic models7. The PCR-L1O method eV. This shows the high predictive power of PCR as the number of
was then benchmarked against experimental and DFT thermo- explicit DFT calculations is reduced 20 times for each new alloy.
chemical data (Fig. 5a)52,53. Our estimates, shown in orange, lie However, the PCA/PCR employed above can be extended to
nicely close to the 1:1 line. The error bars lie between the ±0:25 account for other effects, likely appearing as new components in
eV typical of DFT predictions, as shown by PBE and BEEF-vdW the expansion in Eq. (5). Examples of these effects are coverage
values obtained by other groups (in dark and light gray, contributions54, transition state search18, site coordination22,35,45,
respectively)52,53. and solvation7,55.
Finally, we have used PCR to generate the first full thermo-
chemical database containing SAAs and NSAs. The SAAs were
obtained by replacing one atom by the guest element, whereas the Discussion
NSA14 were generated by substituting either the overlayer or the Statistical learning provides a robust toolbox to go beyond the
sub-surface layer (Fig. 5b–d). The host elements were those listed traditional interpretation of the energetics of intermediates on
in Fig. 4 and the guest elements also included Fe, Co, and Re, for a transition metal surfaces. By applying PCA on a set of formation

NATURE COMMUNICATIONS | (2019)10:4687 | https://doi.org/10.1038/s41467-019-12709-1 | www.nature.com/naturecommunications 5

ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12709-1

energies obtained by DFT, we extracted only two descriptors for 4. Nørskov, J. K., Bligaard, T., Rossmeisl, J. & Christensen, C. H. Towards the
the reactivity of metals and alloys. A close analysis to these computational design of solid catalysts. Nat. Chem. 1, 37–46 (2009).
descriptors revealed their nature as covalent and ionic contribu- 5. Sutton, J. E., Guo, W., Katsoulakis, M. A. & Vlachos, D. G. Effects of
correlated parameters and uncertainty in electronic-structure-based chemical
tions, which can be traced back to the well-known d-band model, kinetic modelling. Nat. Chem. 8, 331–337 (2016).
the reduction potential of the metal, as well as conjugation and 6. Seh, Z. W. et al. Combining theory and experiment in electrocatalysis: insights
conformational effects. Then, PCR enabled us to do a fast survey into materials design. Science 355, eaad4998 (2017).
on the adsorption on SAAs and NSAs, using just three adsorbates 7. Li, Q., García-Muelas, R. & López, N. Microkinetics of alcohol reforming for
as predictors: *O, *OH, and *CCHOH. A minimum of DFT H2 production from a FAIR density functional theory database. Nat.
Commun. 9, 526 (2018).
energy evaluations (around 1800) is required to predict the full set 8. Kitchin, J. R. Machine learning in catalysis. Nat. Catal. 1, 230 (2018).
of 31,000 formation energies with high accuracy, as the error bars 9. Tran, K. & Ulissi, Z. W. Active learning across intermetallics to guide
are comparable to DFT ones. This approach is modular and can discovery of electrocatalysts for CO2 reduction and H2 evolution. Nat. Catal.
be extended to other materials and systems, such as the prediction 1, 696 (2018).
of activation energies for elementary steps or the quantification of 10. Singh, A. K., Montoya, J. H., Gregoire, J. M. & Persson, K. A. Robust and
synthesizable photocatalysts for CO2 reduction: a data-driven materials
solvent and coverage effects, thus paving the way for reliable discovery. Nat. Commun. 10, 443 (2019).
thermochemical models suited to heterogeneous catalysis. 11. Sutton, J. E. & Vlachos, D. G. Building large microkinetic models with first-
principles’ accuracy at reduced computational cost. Chem. Eng. Sci. 121,
Methods 190–199 (2015).
Computational details. We performed DFT calculations with the Vienna Ab- 12. Zaffran, J., Michel, C., Auneau, F., Delbecq, F. & Sautet, P. Linear energy
initio Simulation Package, VASP56,57, for 71 species derived from C1 and C2 relations as predictive tools for polyalcohol catalytic reactivity. ACS Catal. 4,
alcohol decomposition7,49 on closed-packed surfaces of Cu, Ag, Au, Ni, Pd, Pt, Rh, 464–468 (2014).
Ir, Ru, Os, Zn, and Cd. The structures can be retrieved from ioChem-BD58,59, 13. Alonso, D. M., Wettstein, S. G. & Dumesic, J. A. Bimetallic catalysts for
which is a FAIR (Findability, Accessibility, Interoperability, and Reusability) upgrading of biomass to fuels and chemicals. Chem. Soc. Rev. 41, 8075–8098
database43. Instructions about data management are provided in Supplementary (2012).
Methods. The functional of choice was PBE60 with the D2 dispersion corrections of 14. Greeley, J. & Mavrikakis, M. Alloy catalysts designed from first principles. Nat.
Grimme and our reparameterized values for metals61,62. The structural and elec- Mater. 3, 810–815 (2004).
tronic parameters for the metals can be found on Supplementary Tables 1–2. The 15. Nikolla, E., Schwank, J. & Linic, S. Measuring and relating the electronic
present setup follows the current gold standard in DFT calculations42. Core elec- structures of nonmodel supported catalytic materials to their performance. J.
trons were represented by Projector Augmented Wave pseudopotentials63,64 and Am. Chem. Soc. 131, 2747–2754 (2009).
valence electrons were represented by plane waves with a kinetic energy cutoff of 16. Marcinkowski, M. D. et al. Pt/Cu single-atom alloys as coke-resistant catalysts
450 eV. The calculated lattice parameters for the metals show good agreement with for efficient C-H activation. Nat. Chem. 10, 325 (2018).
experimental values, as detailed in the Supplementary Information. Metal surfaces 17. Lucci, F. R. et al. Controlling hydrogen activation, spillover, and desorption
were modeled by four-layer slabs, where the two uppermost layers were fully with Pd-Au single-atom alloys. J. Phys. Chem. Lett. 7, 480–485 (2016).
relaxed and the bottom ones were fixed to the bulk distances. We selected the (111) 18. Li, Z., Wang, S., Chin, W. S., Achenie, L. E. & Xin, H. High-throughput
surfaces for the
pffiffiffi fcc p
metals
ffiffiffi and the (0001) for the hcp ones. The adsorption was screening of bimetallic catalysts enabled by machine learning. J. Mater. Chem.
studied on 2 3 ´ 2 3  R30 supercells. The vacuum between the slabs was set A 5, 24131–24138 (2017).
larger than 13 Å and the dipole correction was applied in z direction65. The 19. Duchesne, P. N. et al. Golden single-atomic-site platinum electrocatalysts.
Brillouin zone was sampled by a Γ-centered 3 ´ 3 ´ 1 k-points mesh generated Nat. Mater. 17, 1033 (2018).
through the Monkhorst–Pack method66. For each species, several conformations 20. Greiner, M. T. et al. Free-atom-like d states in single-atom alloy catalysts. Nat.
were calculated by DFT7,49, but only the most stable ones were taken for sub- Chem. 10, 1008 (2018).
sequent analysis. The gas-phase molecules were relaxed in a cubic box with 20 Å 21. Batchelor, T. A. A. et al. High-entropy alloys as a discovery platform for
sides. For PCA and PCR, the diagonalizations were done with Maple using double electrocatalysis. Joule 3, 834–845 (2019).
precision. 22. Choksi, T. S., Roling, L. T., Streibel, V. & Abild-Pedersen, F. Predicting
adsorption properties of catalytic descriptors on bimetallic nanoalloys with
Data availability site-specific precision. J. Phys. Chem. Lett. 10, 1852–1859 (2019).
All relevant data are available from the authors. The matrices E, E′, and E″ from the 23. Wei, J. N., Duvenaud, D. & Aspuru-Guzik, A. Neural networks for the
training set are uploaded as Supplementary Files matrix-E.csv, matrix-E- prediction of organic chemistry reactions. ACS Cent. Sci. 2, 725–732 (2016).
prime.csv, and matrix-E-second.csv, respectively. The matrix E′ for SAAs and 24. Ulissi, Z. W., Medford, A. J., Bligaard, T. & Nørskov, J. K. To address surface
NSAs is presented on matrix-E-prime-alloys.csv and the 31,453 structures reaction network complexity using scaling relations machine learning and
predicted are listed on alloys-prediction.csv. The matrix files are labeled using DFT calculations. Nat. Commun. 8, 14621 (2017).
a succinct notation detailed in Supplementary Methods and Supplementary Table 5. The 25. O’Connor, N. J., Jonayat, A. S. M., Janik, M. J. & Senftle, T. P. Interaction
structures of all species can be downloaded from ioChem-BD58 following ref. 59. trends between single metal atoms and oxide supports identified with density
functional theory and statistical learning. Nat. Catal. 1, 531 (2018).
26. Butler, K. T., Davies, D. W., Cartwright, H., Isayev, O. & Walsh, A. Machine
Code availability learning for molecular and materials science. Nature 559, 547 (2018).
The LibreOffice Calc spreadsheets, Maple worksheets, and Python scripts are available 27. Chowdhury, A. J. et al. Prediction of adsorption energies for chemical species
from the authors. The data from pure metals and alloys are processed in pca- on metal catalyst surfaces using machine learning. J. Phys. Chem. C 122,
regressions.ods and pca-alloys.ods spreadsheets, respectively. The Maple 28142–28150 (2018).
worksheet used for diagonalizations is diagonalization.mw. The python script that 28. Benson, S. W. Thermochemical Kinetics (Wiley, 1976).
generates all thermochemical data for alloys is pca-alloys.py. 29. Benson, S. W. et al. Additivity rules for the estimation of thermochemical
properties. Chem. Rev. 69, 279–324 (1969).
30. Hammer, B. & Nørskov, J. K. Why gold is the noblest of all the metals. Nature
Received: 25 February 2019; Accepted: 23 August 2019; 376, 238–240 (1995).
31. Hammer, B., Morikawa, Y. & Nørskov, J. K. CO chemisorption at metal
surfaces and overlayers. Phys. Rev. Lett. 76, 2141 (1996).
32. Abild-Pedersen, F. et al. Scaling properties of adsorption energies for
hydrogen-containing molecules on transition-metal surfaces. Phys. Rev. Lett.
99, 016105 (2007).
References 33. Calle-Vallejo, F., Martínez, J. I., García-Lastra, J. M., Rossmeisl, J. & Koper, M.
1. Besson, M., Gallezot, P. & Pinel, C. Conversion of biomass into chemicals over T. M. Physical and chemical nature of the scaling relations between adsorption
metal catalysts. Chem. Rev. 114, 1827–1870 (2013). energies of atoms on metal surfaces. Phys. Rev. Lett. 108, 116103 (2012).
2. Resasco, D. E., Wang, B. & Sabatini, D. Distributed processes for biomass 34. Montemore, M. M. & Medlin, J. W. Site-specific scaling relations for
conversion could aid UN Sustainable Development Goals. Nat. Catal. 1, 731 hydrocarbon adsorption on hexagonal transition metal surfaces. J. Phys.
(2018). Chem. C 117, 20078–20088 (2013).
3. Jones, G. Industrial computational catalysis and its relation to the digital 35. Calle-Vallejo, F. et al. Finding optimal surface sites on heterogeneous catalysts
revolution. Nat. Catal. 1, 311 (2018). by counting nearest neighbors. Science 350, 185–189 (2015).

6 NATURE COMMUNICATIONS | (2019)10:4687 | https://doi.org/10.1038/s41467-019-12709-1 | www.nature.com/naturecommunications

NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12709-1 ARTICLE

36. Stenlid, J. H. & Brinck, T. Extending the σ-hole concept to metals: An 59. García-Muelas, R. Statistical learning goes beyond the d-band model
electrostatic interpretation of the nanostructural effects in gold and platinum providing the thermochemistry of adsorbates on transition metals: Dataset.
catalysis. J. Am. Chem. Soc. 139, 11012–11015 (2017). Stored in ioChem-BD, Ref. [58]. https://doi.org/10.19061/iochem-bd-1-43
37. Salciccioli, M., Chen, Y. & Vlachos, D. G. Density functional theory-derived (2019).
group additivity and linear scaling methods for prediction of oxygenate 60. Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized gradient approximation
stability on metal catalysts: adsorption of open-ring alcohol and polyol made simple. Phys. Rev. Lett. 77, 3865–3868 (1996).
dehydrogenation intermediates on Pt-based metals. J. Phys. Chem. C 114, 61. Grimme, S. Semiempirical GGA-type density functional constructed with a
20155–20166 (2010). long-range dispersion correction. J. Comput. Chem. 27, 1787–1799 (2006).
38. Salciccioli, M., Edie, S. M. & Vlachos, D. G. Adsorption of acid, ester, and 62. Almora-Barrios, N., Carchini, G., Błoński, P. & López, N. Costless derivation
ether functional groups on Pt: fast prediction of thermochemical properties of of dispersion coefficients for metal surfaces. J. Chem. Theory Comput. 10,
adsorbed oxygenates via DFT-based group additivity methods. J. Phys. Chem. 5002–5009 (2014).
C 116, 1873–1886 (2012). 63. Blöchl, P. E. Projector augmented-wave method. Phys. Rev. B 50, 17953–17979
39. Suntivich, J. et al. Design principles for oxygen-reduction activity on (1994).
perovskite oxide catalysts for fuel cells and metal-air batteries. Nat. Chem. 3, 64. Kresse, G. & Joubert, D. From ultrasoft pseudopotentials to the projector
546–550 (2011). augmented-wave method. Phys. Rev. B 59, 1758–1775 (1999).
40. Suntivich, J., May, K. J., Gasteiger, H. A., Goodenough, J. B. & Shao-Horn, Y. 65. Makov, G. & Payne, M. C. Periodic boundary conditions in ab initio
A perovskite oxide optimized for oxygen evolution catalysis from molecular calculations. Phys. Rev. B 51, 4014–4022 (1995).
orbital principles. Science 334, 1383–1385 (2011). 66. Monkhorst, H. J. & Pack, J. D. Special points for Brillouin-zone integrations.
41. Hong, W. T. et al. Toward the rational design of non-precious transition metal Phys. Rev. B 13, 5188–5192 (1976).
oxides for oxygen electrocatalysis. Energy Environ. Sci. 8, 1404–1427 (2015).
42. Lejaeghere, K. et al. Reproducibility in density functional theory calculations
of solids. Science 351, 1415 (2016). Acknowledgements
43. Wilkinson, M. D. et al. The FAIR guiding principles for scientific data We thank AGAUR 2017 SGR 90 and MCIU RTI2018-101394-B-I00 projects for financial
management and stewardship. Sci. Data 3, 160018 (2016). support. The Barcelona Supercomputing Center – MareNostrum (BSC-RES) is
44. Yao, K., Herr, J. E., Brown, S. N. & Parkhill, J. Intrinsic bond energies from a acknowledged for providing generous computer resources. We thank Professors C. Bo,
bonds-in-molecules neural network. J. Phys. Chem. Lett. 8, 2689–2694 (2017). R. Guimerà, M. Sales-Pardo, F. Abild-Pedersen, and N. Almora-Barrios for fruitful
45. Andersen, M., Levchenko, S., Scheffler, M. & Reuter, K. Beyond scaling discussions.
relations for the description of catalytic materials. ACS Catal. 9, 2752–2759
(2019). Author contributions
46. Xin, H. & Linic, S. Communications: exceptions to the d-band model of R.G.-M. performed the numerical calculations. R.G.-M. and N.L. contributed to analyze
chemisorption on metal surfaces: the dominant role of repulsion between the data and to write the manuscript.
adsorbate states and metal d-states. J. Chem. Phys. 132, 221101 (2010).
47. Janet, J. P. & Kulik, H. J. Resolving transition metal chemical space: Feature
selection for machine learning and structure-property relationships. J. Phys. Competing interests
Chem. A 121, 8939–8954 (2017). The authors declare no competing interests.
48. Friedman, J., Hastie, T. & Tibshirani, R. in The elements of statistical learning,
2nd edn Pages: xii; 79–80; 408–409; 534–541 (Springer Series in Statistics,
2001).
Additional information
Supplementary information is available for this paper at https://doi.org/10.1038/s41467-
49. García-Muelas, R., Li, Q. & López, N. Density functional theory comparison of
019-12709-1.
methanol decomposition and reverse reactions on metal surfaces. ACS Catal.
5, 1027–1036 (2015).
Correspondence and requests for materials should be addressed to R.G.-M. or N.L.
50. Lide, D. CRC Handbook of Chemistry and Physics 84th edn (CRC Press LLC,
2003–2004).
Peer review information Nature Communications thanks the anonymous reviewer(s) for
51. Montemore, M. M. & Medlin, J. W. A unified picture of adsorption on
their contribution to the peer review of this work. Peer reviewer reports are available.
transition metals through different atoms. J. Am. Chem. Soc. 136, 9272–9275
(2014).
Reprints and permission information is available at http://www.nature.com/reprints
52. Wellendorff, J. et al. A benchmark database for adsorption bond energies to
transition metal surfaces and comparison to selected DFT functionals. Surf.
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in
Sci. 640, 36–44 (2015).
published maps and institutional affiliations.
53. Silbaugh, T. L. & Campbell, C. T. Energies of formation reactions measured
for adsorbates on late transition metal surfaces. J. Phys. Chem. C 120,
25161–25172 (2016).
54. Xu, Z. & Kitchin, J. R. Probing the coverage dependence of site and adsorbate Open Access This article is licensed under a Creative Commons
configurational correlations on (111) surfaces of late transition metals. J. Phys. Attribution 4.0 International License, which permits use, sharing,
Chem. C 118, 25597–25602 (2014). adaptation, distribution and reproduction in any medium or format, as long as you give
55. Garcia-Ratés, M., García-Muelas, R. & López, N. Solvation effects on appropriate credit to the original author(s) and the source, provide a link to the Creative
methanol decomposition on Pd(111), Pt(111), and Ru(0001). J. Phys. Chem. C Commons license, and indicate if changes were made. The images or other third party
121, 13803–13809 (2017). material in this article are included in the article’s Creative Commons license, unless
56. Kresse, G. & Furthmüller, J. Efficiency of ab-initio total energy calculations for indicated otherwise in a credit line to the material. If material is not included in the
metals and semiconductors using a plane-wave basis set. Comput. Mater. Sci. article’s Creative Commons license and your intended use is not permitted by statutory
6, 15–50 (1996). regulation or exceeds the permitted use, you will need to obtain permission directly from
57. Kresse, G. & Furthmüller, J. Efficient iterative schemes for ab initio total- the copyright holder. To view a copy of this license, visit http://creativecommons.org/
energy calculations using a plane-wave basis set. Phys. Rev. B 54, 11169–11186 licenses/by/4.0/.
(1996).
58. ÁlvarezMoreno, M. et al. Managing the computational chemistry big data
problem: The ioChem-BD platform. J. Chem. Inf. Model. 55, 95–103 (2015). © The Author(s) 2019

NATURE COMMUNICATIONS | (2019)10:4687 | https://doi.org/10.1038/s41467-019-12709-1 | www.nature.com/naturecommunications 7