Role of monodentate formate in product selectivity for CO₂ hydrogenation on Pd-based alloy catalysts

Abstract

The hydrogenation of CO₂ to methanol is a key reaction for sustainable fuel synthesis. A crucial aspect of the catalytic mechanism is the role of monodentate formate (HCOO^m*) in the initial steps of CO₂ hydrogenation on Pd-based alloy catalysts, which we have investigated using density functional theory (DFT) together with subgroup discovery (SGD) analysis. The reactivity and stability of CO₂ and formate species are examined on monometallic Pd, Cu, Zn surfaces and alloyed CuPd and PdZn surfaces. PdZn surfaces show low activation energy barriers for CO₂ hydrogenation and, combined with weak CO₂^δ− adsorption energy, this suggests that an Eley–Rideal mechanism may dominate over Langmuir–Hinshelwood pathways. The adsorption energy of the monodentate formate intermediate is found to correlate significantly with the activation energy of CO₂ hydrogenation across all investigated facets, where stronger adsorption yields lower activation energy, enabling its use as a predictive descriptor. To determine possible new catalytic materials, a dataset of 49 Pd-based single-atom alloy (SAA) surfaces is screened using SGD, identifying key electronic parameters, the dopant and site electron affinity, as drivers of strong HCOO^m* adsorption. The obtained subgroup discovery rules highlight Mo, Nb, and W as promising earth-abundant dopants for Pd-based catalysts, further confirmed by DFT calculations. These findings offer a mechanistic rationale for catalyst design and demonstrate the utility of AI-guided screening in identifying efficient, sustainable alloy compositions to be used as catalysts for CO₂ conversion.

---

1 Introduction

The hydrogenation of CO₂ to produce methanol is a promising approach to mitigate industrial emissions of a greenhouse gas while simultaneously synthesising “green” fuel and a chemical substrate of considerable societal value. The commercially-leading catalyst for hydrogenation of CO₂ is a Cu/ZnO/Al₂O₃ (CZA) material, with both the ZnO support and the metallic Cu playing key roles in the catalytic mechanism.¹ CZA is a very efficient catalyst but must operate under high temperature and pressure conditions, and is prone to deactivation in the presence of the water produced in the reaction.² Recently, Pd-based supported catalysts have been investigated for selective hydrogenation of CO₂ to methanol, showing significantly improved water tolerance compared to the prior catalysts and desirable selectivity to only C₁ products,^3,4 though with undesirable methane produced as well as the desired methanol product.

To explain the reactivity of Cu and Pd catalysts for CO₂ hydrogenation, experimental studies have been beneficially complemented by computational simulations using density functional theory (DFT).^5–9 The commonly modelled low energy (111) surfaces of face-centred-cubic (FCC) Cu and Pd catalysts show high energy barriers in reactions involving carbon dioxide and formate, an important intermediate in methanol synthesis, with less stable (100) and (110) facets likely to be more reactive and catalytically relevant.^7,8,10 On Cu and Pd surfaces, CO₂ has been shown to adsorb in either physisorbed (yielding CO₂) or chemisorbed (yielding CO₂^δ−) configurations.^7,8 Physisorbed CO₂ orients parallel to the surface at a distance of over 3 Å and no chemical bond forms between the catalyst and the adsorbate, which means the interaction is only marginally exothermic.^7,8 In contrast, the chemisorbed CO₂^δ− assumes a bidentate geometry with the formation of C–metal and O–metal bonds, and electron density transfers from the metal surface to the π* orbital of the carbon atom, resulting in a distorted geometry with a reduction of the O–C–O bond angle (to 140°) that helps to stabilise the negative charge (^δ−).^7,8 On Zn (0001) surfaces, only the weakly bonded physisorbed CO₂ species forms, with chemisorption on Zn surfaces uncommon in the absence of surface oxygen or a high electrostatic potential.^11–15

The activity of Pd surfaces for CO₂ hydrogenation is limited by the highly endothermic steps in the formation of the formate intermediate, with subsequent hydrogenation of formate to either methane or methanol having lower activation energies.^8,9 The initial stages of electrocatalytic CO₂ hydrogenation on Pd-based catalysts have been shown to involve a monodentate formate (HCOO^m) intermediate, with only a single metal–oxygen bond between the catalyst and the molecule.¹⁶ The so-called formate pathway, where the carbon in CO₂ is hydrogenated to yield the HCOO intermediate, has been linked to an increased selectivity towards methanol in catalytic systems that involve Cu, Pd and Zn metals.^7,17–19

Alloying Pd catalysts with Zn has been shown to facilitate tuning of the reaction selectivity, with no methane observed in thermal CO₂ hydrogenation at ambient pressure and temperatures up to 500 K.^20,21 X-ray and neutron powder diffraction measurements were used to describe precisely the structure of the PdZn catalyst, identifying prominence of a tetragonal β₁-PdZn alloy phase.²² The structure of the PdZn nanoparticle catalysts was studied using various microscopy methods and identified tetragonal crystals with greatest abundance of (101) and (110) facets.^23–26 The 1 : 1 PdZn alloy was identified, showing features of a perfectly ordered material with a pattern of alternating stripes of Pd and Zn when viewed in the xy-plane.^23–27 Theoretical studies of the PdZn alloy at varying composition confirmed the 1 : 1 ratio as most energetically favourable, and identified the most stable structure as body-centred tetragonal, in agreement with experiment.²⁸

An alternative catalyst composition of Pd alloyed with Cu has been shown to also enhance CO₂ hydrogenation reactivity, relative to monometallic systems, with geometry²⁹ and electronic³⁰ effects from the CuPd alloying important in defining the selectivity towards C₁ and C₂ products. Alloying Pd with Cu enhances the thermally-driven CO₂ hydrogenation reaction, by means of a synergetic effect of hydrogen spill-over from Pd onto the Cu active sites,³¹ as corroborated by multiple studies.^{27,30,32–34} Challenges remain, however, with respect to understanding the reaction mechanism on these surfaces, particularly with respect to the stability of the formate species and the most prominent reaction mechanism.³⁵

The mechanism of CO₂ hydrogenation to methanol on these supported alloy nanoparticle catalysts is demonstrably linked to the formate intermediate,³⁵ and application of novel computational tools to probe the chemical reactivity of this key intermediate can facilitate new understanding and rational catalyst design. Simulations of homogeneous alloys can be coupled with dilute-limit representations of the catalyst, i.e., single-atom alloy (SAA) catalysts, to understand concentration effects; furthermore, SAAs offer a platform for high-throughput screening of catalysts, for tuning of the active sites, exploring structure–activity relationships, and improving reactivity with small metal loadings.^34,36–38 Theoretical approaches have successfully explained the chemical phenomena occurring on SAA catalysts, such as enhancement or creation of new active sites.^37–41 Additionally, the surface structure of SAA catalysts is typically constrained to that of the host material, making the modelling of SAA surfaces more accessible than more disordered alloy systems, as the number of surface configurations are greatly reduced; the consequence is more efficient theoretical probing of catalytic materials can be performed on a quantitative basis.³⁸ Thus, SAA catalysts offer significant opportunities to strengthen understanding of reaction mechanisms, and to tailor active sites to improve efficacy for CO₂ hydrogenation.⁴²

Here, we present an extensive study of the stability and reactivity of CO₂ and formate species over important surfaces for Pd-based catalysts with the aim of elucidating reaction details and progressing catalyst design. We apply density functional theory (DFT) to calculate the stability and reactivity of CO₂ and hydrogen on homogeneously alloyed body-centred tetragonal (BCT) PdZn (101) and (110) surfaces, and a body-centred cubic (BCC) CuPd (110) surface. The hydrogenation of CO₂ to produce formate is explicitly considered, and the performance of alloys is referenced against their monometallic counterparts (Pd, Zn, and Cu). We build on our mechanistic results by surveying the stability of formate on SAA surfaces of Pd (100), (110), (111), and (211) doped with Co, Cu, Ga, Ir, Ni, Os, Pt, Rh, Ru or Zn. The SAA surfaces have been described using 13 descriptive parameters intrinsic to the materials, resulting in a complex dataset that warrants the use of Artificial Intelligence (AI) to identify the most chemically relevant information. The method applied herein is subgroup discovery (SGD), a focused AI approach devoted to finding subsets or subgroups displaying exceptional patterns within a data set.⁴³ SGD is a powerful and versatile tool that can be applied to both scientific and non-scientific problems where the target property of interest is well defined either numerically or categorically. Through the combined use of DFT and SGD, we are able to highlight and confirm novel and earth-abundant transition metal alloy candidates with potential to increase the efficacy of Pd-based catalysts for CO₂ hydrogenation.

2 Computational methods

2.1 Electronic structure calculations

The Fritz Haber Institute ab initio materials simulations (FHI-aims) software package has been used for full potential all-electron DFT calculations, with the Pythonic Atomic Simulation Environment (ASE) used for management of calculation geometries.^44,45 The default convergence criteria within FHI-aims for self-consistent field (SCF) calculations were used, i.e., the changes between the current and previous SCF iterations in charge density, sum of eigenvalues, and total energy were below N × 1.67 × 10⁻⁵ e a₀⁻³, 10⁻³ eV, and 10⁻⁶ eV, respectively, where N is the number of atoms in the model. Scalar relativistic treatment of kinetic energy for all elements was achieved by the atomic zero-order regular approximation (ZORA),⁴⁶ and a Gaussian-type broadening with width of 0.01 eV was applied to the occupation of electronic states. The accurate semi-local Bayesian error estimation density functional (mBEEF)⁴⁷ exchange–correlation density functional (DF), as implemented in the LibXC DF library,⁴⁸ has been previously identified as a suitable choice for modelling alloy Cu, Pd and Zn surfaces,¹⁰ and this DF has been applied throughout. A reciprocal space sampling of one k-point per (0.018 × 2π) Å⁻¹, and a default “light” basis set (version: 2020), have been used for geometry optimisations, providing structural accuracy and electronic accuracy in Mulliken charge population analysis.^10,44,49 Spin-paired calculations were used in periodic calculations; gas-phase adsorbate structures were calculated with both spin-paired and spin-unpaired configurations, and the energy of the more stable system was considered for reference to periodic calculations. For geometry optimisations, convergence was deemed complete when forces on all unconstrained atoms were less than 0.01 eV Å⁻¹; a force correction term was applied for geometry optimisations also, which suppresses the strong fluctuations in the Hellmann–Feynman forces when the electron densities are loosely converged and results in speed up of ∼30%, without compromising calculation accuracy.⁴²

Monometallic surfaces were constructed from face-centred cubic (FCC) bulk unit cells of Cu and Pd, and a hexagonal close-packed (HCP) unit cell of Zn, based on the outcomes of the work by Kabalan et al.¹⁰ Supercells of the surfaces were created for Pd (111), (100), and (110); Cu (111), (100), and (110); and Zn (0001), with dimensions of (3 × 3 × 7) unit cells and a vacuum layer of 40 Å added in the z-direction. The bottom three atomic layers were constrained to their bulk positions, with the 4 top surface layers remaining unconstrained; and, due to the one-sided nature of the slab models, a dipole-correction was used in all calculations. The alloy surfaces were prepared similarly to the monometallic surfaces; however, as BCT and BCC primitive cells contained two atoms, (2 × 3 × 7) supercells were used. For the body-centred tetragonal (BCT) PdZn (101) and (110) surfaces, the bimetallic alloy ordering was matched to experimental studies by Peterson et al. and Nowicka et al.,^22,26 with the greater stability of PdZn (110) confirmed by our calculations. The ordering of the body-centred cubic (BCC) CuPd (110) surface was based on a study by Yamauchi et al.⁵⁰ Following the naming convention of Crawley et al. the investigated ordered (o-) alloy catalysts could be described as body-centred tetragonal (BCT) o-Pd_0.5Zn_0.5 and body-centred cubic (BCC) o-Cu_0.5Pd_0.5.⁵¹ Here, we refer to these alloys as PdZn and CuPd, respectively, with crystal structure, composition and ordering implied for simplicity. All prepared facets are presented in Fig. S1 in the SI.

To study the catalysed surface reactions, the adsorption energy (E_ads) measuring the interaction between a surface and adsorbate is deduced from comparison of energies of optimised gas-phase adsorbate (E_adsorbate), pristine surface (E_pristine), and the optimised combined system (E_slab):

E_ads = E_slab − (E_adsorbate + E_pristine)

For the formate species, the reference E_adsorbate is instead taken from ½H₂ and CO₂ to avoid referencing against the unstable gas-phase formate radical. The definition used is such that negative E_ads indicates an exothermic interaction. The E_ads of formate species on the catalyst surface (*) is considered in monodentate (HCOO^m*) or bidentate (HCOO^b*) form, where ^m and ^b superscripts denote mono- or bi-dentate configurations; carbon dioxide is considered either as physisorbed (CO₂*) or chemisorbed (CO₂^δ−*); and hydrogen as dissociatively adsorbed only (H*), where * indicates a surface-bound species in all cases. To rationalise how the alloying of Pd, Cu, and Zn changes the reaction profile for CO₂ hydrogenation to formate, models were constructed for the hydrogenation of physisorbed and chemisorbed CO₂ to HCOO^m* and HCOO^b* species. The lowest lying transition state (TS) is reported in each case, as described in more detail in Section S3 of the SI. The final product state taken as an energy reference corresponds to the strongly surface-bound HCOO^b* intermediate. For the calculation of kinetic barriers, the climbing-image MLNEB method has been used to identify saddle points and minimum energy paths (MEP).^52,53 In line with our previous work, a spring constant of 0.05 eV Å⁻¹ has been used throughout; the convergence criteria is for forces on all unconstrained atoms to be below 0.05 eV Å⁻¹, with energy uncertainty below 0.03 eV.⁸

For the single atom alloy (SAA) catalysts, calculations of Pd (111), (100), (110), and (211) (Fig. S2 in the SI) surfaces were considered with a single atom in the top atomic layer substituted by one of the following elements: Co, Cu, Ga, Ir, Ni, Os, Pt, Rh, Ru and Zn. The structures were all subjected to geometry relaxation,⁵⁴ with the self-consistent field (SCF) parameters consistent with the above details except for structures containing Co and Ni, where a spin-unpaired calculation was chosen, with the default initial magnetic moment set to 0.01 so as to allow dynamic spin convergence; for these cases, the calculations also required the occupation broadening was increased 0.1 eV, the mixer was explicitly set to Pulay, and the Kerker preconditioner was disabled, to aid SCF convergence. For studies of the SAA catalysts, the initial structure of monodentate formate species were adapted from chemisorbed CO₂ geometries on the Pd SAA facets,⁴² and refined with the MACE-MP forcefield with a pre-trained “small” model and a van der Waals correction.^55–57 Full details of the calculation settings are provided in the SI.

2.2 Subgroup discovery

From the structures optimised with DFT, a set of basic candidate parameters characterising the (clean) adsorption sites and the formate adsorption energy values were obtained. This data served as inputs for the SGD approach. The dataset constitutes 49 symmetrically inequivalent monodentate formate structures adsorbed on the surfaces of SAAs. The goal of the SGD analysis is to identify the key parameters that describe the adsorption properties of monodentate formate. A detailed description of SGD can be found in the work of S. Wrobel,⁴³ with contextual examples provided by B. R. Goldsmith et al.⁵⁸ and Foppa et al.⁵⁹ Briefly, the SGD algorithm identifies descriptions of subsets of data (subgroups) displaying outstanding distributions of a given target of interest, in our case the adsorption energy of HCOO^m*. The SGD descriptions or rules are typically inequalities constraining the value of key basic parameters, selected from all the offered candidate parameters. The rules provided by SGD can be used to identify new surface sites with the desired adsorption properties. The SG identification is based on the maximization of a quality function Q(S, P), which ensures a balance between the size of the SG and its utility, e.g., the preference for low or high values of the target of interest. Here, we utilize the negative mean shift to identify SGs of surface sites with low adsorption energy of HCOO^m*.

The 13 candidate descriptive parameters considered for SGD in this work correspond to two classes: SA and surface adsorption site. These parameters are chosen to comprise the SAA surfaces’ electronic and geometric properties. The SA parameters are the free-atom Pauling electronegativity (PE_SA), the electron affinity (EA_SA), the ionization potential (IP_SA), and the radii of the s-, p-, d-, and valence orbitals (r_s-SA, r_p-SA, r_d-SA, r_val-SA). Surface sites are described using PE_site, IP_site, and EA_site, where the subscript “site” implies the average of these properties over all the atoms in the adsorption ensemble, i.e. nearest neighbouring metal atoms. The geometry of the adsorption site is included through the coordination number (CN), the generalized CN (gen-CN), and the number of atoms in the adsorption site (N_site). The technical details for the SGD studies are provided in the Section S5 in the SI.

3 Results and discussion

3.1 CO₂ adsorption on metallic and 1 : 1 alloy surfaces

To gain insight into the interaction of the key CO₂ hydrogenation reaction substrates with the metallic Cu, Pd, Zn and CuPd and PdZn alloy surfaces, CO₂ adsorption and co-adsorption with H was investigated on (111), (100), (110) surfaces of Cu and Pd, the (0001) Zn facet, CuPd (110), and (110) and (101) surfaces of PdZn. The most stable calculated adsorption energy, E_ads, for CO₂ on monometallic and alloy surfaces are given in Fig. 1. The physisorbed CO₂ species has very similar stability on all investigated surfaces, with E_ads of between −0.17 and −0.10 eV; thus, the reactivity of CO₂ on specific surfaces is difficult to deduce based on these interactions. In contrast, the stability of chemisorbed CO₂^δ− is more varied and carries more information about reactivity.

Fig. 1 E_ads of CO₂ and CO₂^δ− in units of eV, with and without co-adsorbed H in the most favourable adsorption position, across Pd (111), (100), (110), Zn (0001), Cu (111), (100), (110), CuPd (110) and PdZn (101) and (110) surfaces. No metastable structures of CO₂^δ− were found on Cu (111) and Zn (0001) surfaces. The numerical data is available in Table S4 in the SI.

The exothermic CO₂^δ− adsorption on the Pd (100) and (110) surfaces implies that these surfaces may chemically activate CO₂, with the latter shown in Fig. 2, but the interaction is weak and therefore any chemisorbed state may be short-lived, matching our previous conclusions obtained using a PBE+TS method.⁸ The E_ads of chemisorbed CO₂^δ− on the denser Cu and Pd surfaces was calculated as 0.39 eV and 0.15 eV on Cu (100) and Pd (111), respectively, and no CO₂^δ− structure was observed on Cu (111), in line with observations of Higham et al.⁷ The calculated E_ads of CO₂^δ− on Cu (110) in the current work is weaker than reported by Higham et al., which is attributed to differences in calculations settings; we note that the mBEEF exchange–correlation density functional used here has proven greater accuracy for both material and adsorption processes in recent testing and other works.^8,10,60

Fig. 2 The CO₂ adsorption geometries on Pd (110), with subfigures 1, and 2 corresponding to the CO₂ physisorbed along the short-bridge site and CO₂^δ− along the long-bridge site, respectively. The view is along the xy-, xz-, and xz-axis in sub figures A, B, and C, respectively. Subfigures 1C and 2C apply a ball and stick model to highlight bonding interactions. Subfigure 3 presents a magnified and annotate 2C figure. The black lines represent the boundaries of the unit cell. The grey, red, and blue spheres represent the C, O, and Pd atoms, respectively. The numerical data is available in Table S4 in the SI.

The E_ads of CO₂^δ− on PdZn (101) and (110) is 0.22 eV and 0.18 eV, respectively, which is weaker than the strongest adsorption on pure Pd (111). The most stable CO₂^δ− geometries are on the Cu (110), Pd (110), CuPd (110), and PdZn (110) facets for the different metal surfaces, with an E_ads of 0.30 eV, −0.25 eV, 0.18 eV and 0.09 eV, respectively. The adsorbed structures are shown in Fig. 2 and 3, with additional structural information available in Table S4 in the SI.

Fig. 3 The CO₂^δ− adsorption geometries on (1) Cu (110), (2) CuPd (110) and (3) PdZn (110). The view is along the xy-, yz-, and xz-axis in A, B, and C, respectively. The black lines represent the boundaries of the unit cell. The grey, brown, red, blue, and purple spheres represent the C, Cu, O, Pd, and Zn atoms, respectively.

The endothermic CO₂^δ− adsorption on Cu (100), (110), Pd (111), CuPd (110) and PdZn (101) and (110) surfaces shows that CO₂^δ− on these surfaces is likely a transient species. For the CuPd and PdZn alloy surfaces, only a single chemisorbed structure of CO₂^δ− is identified. On the bimetallic PdZn surfaces, the bidentate CO₂^δ− species binds on intermetallic bridging sites, with C more strongly bound at Pd and O at Zn; the preferable geometry of CO₂^δ− adsorbed on CuPd (110) is similar, with Pd–C and Cu–O bonds, and the monometallic adsorption sites are unstable.

Hydrogenation of adsorbed CO₂ to produce formate requires a H* species in the vicinity of the carbon dioxide; therefore, the relative stability of the CO₂ and CO₂^δ− were assessed with co-adsorption of hydrogen, and the results are also provided in Fig. 1. Other than for the Pd surfaces, the initial co-adsorbed CO₂^δ−* and H* state, required for the Langmuir–Hinshelwood mechanism, is exothermic only on CuPd (110), with all other surfaces endothermic or unstable [Zn (0001); Cu (111)]. Given the endothermicity of the Langmuir–Hinshelwood reactant state, alternative mechanisms, such as the Eley–Rideal, might be playing a significant role in initial CO₂ hydrogenation to formate on PdZn surfaces, as was previously observed for Cu surfaces.⁶¹

3.2 HCOO

Three binding modes of the HCOO* intermediate have previously been identified on Cu surfaces, which correspond to two monodentate configurations (HCOO^m*, with and without hydrogen atom interacting with the surface) and a bidentate configuration (HCOO^b*).^62–66 As the FCC Cu and Pd (110) surfaces were previously highlighted as more catalytically active towards CO₂ hydrogenation than the more stable (111) and (100) facets,^7,8 the relevance of the formate species in CO₂ hydrogenation was investigated by studying its stability on the Cu and Pd (110) surfaces, as well as the prominent alloy surfaces of CuPd⁵⁰ [(110)] and PdZn²⁶ [(101) and (110)]. The formate adsorbate was initially positioned 2 Å above the adsorption sites, with either both oxygen atoms, one oxygen and one hydrogen atom, or just one oxygen atom facing towards the surface, with example resulting structures on the PdZn (101) surface depicted in Fig. 4. The E_ads of the HCOO* adsorption modes on Cu (110), Pd (110), CuPd (110), PdZn (101) and PdZn (110) surfaces are shown in Table 1.

Fig. 4 Three types of formate intermediates considered on the BCT PdZn (110) facet: bidentate (1) and monodentate with (2) and without (3) a significant hydrogen–metal interaction, in xy- and yz-view (A and B, respectively); blue, purple, red, grey and white spheres represent Pd, Zn, O, C, and H atoms, respectively; the black lines are boundaries of the simulation cell. The numerical data is available in Table S4.

Table 1 E_ads(HCOO) on FCC Cu and Pd (110) surfaces, BCC CuPd (110) and BCT PdZn (101) and (110) surfaces, relative to gas-phase reactant species (½H₂ and CO₂); ^m* corresponds to a monodentate adsorption geometry and ^b* corresponds to a bidentate adsorption geometry. In some cases, metastable formate configurations were identified with the H moiety directed towards the surface to form a H–M interaction (where M is a metal species), and the E_ads in these cases is given in bold. Models showing unstable adsorption are marked with dashes

	Eads(HCOO^m*)/eV		Eads(HCOO^b*)/eV
Adsorption site	Cu (110)	Pd (110)	Cu (110)	Pd (110)
Short-bridge	−0.40, 0.09	0.14, 0.17	−1.08	−0.81
Long-bridge	—	—	−1.02	−0.65
Adsorption site	CuPd (110)
Pd–Pd bridge	0.58		−0.34
Pd–Cu bridge	0.50		−0.43
Cu–Cu bridge	0.15		−0.50
Adsorption site	PdZn (101)
Pd–Zn–Zn	−0.21		−0.70
Pd–Pd–Zn	−0.28, 0.12		−0.64
Pd–Zn	—		−0.82
Adsorption site	PdZn (110)
Pd–Zn–Zn	0.17		—
Zn–Zn	—		−0.72, −0.75
Pd–Zn	−0.37, −0.03		−0.87
Pd–Pd	—		−0.41, −0.44

---

Notably, the adsorption of monodentate HCOO (HCOO^m*) without a H–M interaction (where M is a metal species) is endothermic on FCC Pd (110) and BCC CuPd (110) surfaces, and exothermic on the FCC Cu (110), BCT PdZn (101) and (110) surfaces. The only exothermic HCOO^m* species with a H–M interaction, which might be produced by initial CO₂ hydrogenation, is for the PdZn (110) surface. The bidentate geometry (HCOO^b*) is more stable for all systems, especially on the short-bridge site of the FCC (110) surfaces of Cu and Pd, with E_ads of −1.08 eV and −0.81 eV, respectively; HCOO^b* is most stable at the Cu–Cu bridge adsorption site of BCC CuPd (110) surface, with E_ads of −0.50 eV. HCOO^b* is more stable on Cu (110) than on Pd (110), and the preferred adsorption on the Cu atoms of BCC CuPd (110) surface layer matches this observation: E_ads of HCOO^b* is −0.34 eV when both oxygens are bound to Pd, and −0.43 eV and −0.50 when bound to one and two Cu atoms, respectively. In contrast, adsorption of HCOO^b* involving multiple Pd or multiple Zn atoms on PdZn (101) and (110) facets is weaker than on the most stable bimetallic Pd–Zn bridge sites, with E_ads of −0.82 and −0.87 eV on PdZn (101) and (110), respectively, i.e., similar to E_ads of HCOO^b* on the Pd (110) surface. The most stable HCOO^b* structure on the BCC CuPd (110) surface is significantly less stable than on the monometallic Cu (110) and Pd (110) surfaces, by 0.58 eV and 0.31 eV, respectively, suggesting that the CuPd catalyst surface is less prone than Cu to poisoning with the formate species.

Under reaction conditions, the formate species is likely to primarily bond in the bidentate configuration, due to its greater stability on each investigated metal and alloy surface. However, when a HCOO species is formed from CO₂ and H*, the H* transfer from the surface means a non-vertical or non-bidentate position must initially form to facilitate the transfer of the surface-bound hydrogen onto the carbon; this reaction step is highly endothermic, as previously shown on Cu and Pd facets.^7,8,65 Horizontally positioned formate species that correspond to TS structures observed in simulations on Cu⁷ and Pd⁸ surfaces are not metastable; however, the observed metastable HCOO^m* species with an H–M interaction resembles the TS structure for HCOO formation in an alternative study of formate desorption on Cu (110), and is proposed to be an important intermediate.⁶⁵

Mulliken charge analysis of the HCOO configurations shows that a major part of the electron density withdrawn from the surfaces is towards the oxygen atoms. The C atom has a charge that is 0.20 e higher in the HCOO^b* and HCOO^m* formate structure, with H–M interaction, than for the physisorbed CO₂. The charge on the H atom in the H–M configuration of HCOO^m* is −0.04 to −0.08 e lower (i.e., less positive charge, more electrons) than in HCOO^b*. The H–M interaction may enhance an Eley–Rideal-type mechanism for the formation of formate, where physisorbed CO₂ interacts with the surface-adsorbed H, thus bypassing the CO₂ chemisorption necessary for the Langmuir–Hinshelwood mechanism.

The stability of the formate adsorbates was explored further, considering the most stable HCOO^b* and the HCOO^m* intermediates with an H–M bond that might enable alternative reaction mechanisms. The HCOO^b* and HCOO^m* intermediates were modelled on the FCC Pd (111), (100), FCC Cu (111), (100), HCP Zn (0001), and the results compare to those presented in Table 1. The E_ads of the HCOO^b* and HCOO^m* with a H–M interaction were calculated on each surface and are shown in Fig. 5. No metastable HCOO^m* structure has been identified on the Zn (0001) surface, similar to the inability to chemisorb CO₂^δ−.

Fig. 5 E_ads of surface-bound (*) monodentate (^m) and bidentate (^b) formate (HCOO) intermediates on FCC Pd (111), (100), (110), Cu (111), (100), (110), HCP Zn (0001), BCC CuPd (110), and PdZn (101) and (110) surfaces. HCOO^b* is in the most stable conformation with both oxygen atoms bound to metal atoms, while the HCOO^m* is bound to the surface via a single oxygen atom and a significant H–metal electronic interaction. No metastable HCOO^m* structure has been identified on the Zn (0001) surface. The numerical data is available in Table S4 in the SI.

Li et al. previously reported the adsorption energy of monodentate formate on Pd (111) and Pd (100) surfaces to be 0.76 eV and 0.78 eV weaker than the bidentate species,⁶⁷ respectively, which is comparable to the energy differences of 0.91 eV calculated for both surfaces here. The only exothermic E_ads(HCOO^m*) with a H–M interaction is for the PdZn (110) surface, with a marginally negative E_ads of −0.03 eV. The recognition of the E_ads(HCOO^m*) outlier informs both further reactivity investigations and hints at potential targets for SGD analysis.

3.3 Reaction energetics

The mechanism of CO₂ hydrogenation on Cu (100), (110) and (111) facets has been proposed to involve the monodentate HCOO^m* species, from which the HCOO^b* is formed without an energy barrier.⁶⁶ Notably, the calculated E_ads of the TS for the initial CO₂ hydrogenation to formate is similar to the adsorption energy of the HCOO^m* species.⁶⁶ An alternative mechanism, with a horizontally positioned HCOO^b* species, including the H–M interaction and two O–M bonds, was recently reported from DFT calculations with a barrierless rotation to HCOO^b* on Cu (100) and (110) facets.⁷ The HCOO^m* species has been highlighted as a key intermediate in electrochemical reduction of CO₂ on Pd surfaces, and so there is a clear need to investigate further.¹⁶ The activation energy of hydrogenation of physisorbed CO₂ to formate, E_a(CO₂ + H), activation energy of CO₂^δ− hydrogenation to formate, E_a(CO₂^δ− + H), and the reverse reaction, i.e., the activation energy for decomposition of bidentate formate, E_a(HCOO^b) were calculated, and the reaction energy profile is presented in Fig. 6, and the geometric data in Table S4 in the SI.

Fig. 6 Reaction energy profiles of initial hydrogenation of CO₂ to HCOO on FCC Cu, Pd, HCP Zn, BCT PdZn and BCC CuPd surfaces in units of eV. The preferred adsorbed CO₂ geometry for the transition state calculation was the bent chemisorbed CO₂^δ− (solid line). Metastable CO₂^δ− geometry was not found on FCC Cu (111) and HCP Zn (0001), thus the physisorbed CO₂ was used instead (dotted lines); * denotes a surface and surface-bound intermediates. The numerical data is available in Table S4 in the SI.

Considering the TS structures obtained across the investigated Pd, Zn and PdZn surfaces, the carbon–metal bond (d_TS(C–M)) is longer on PdZn than on Pd and the C–H distance in the TS changes significantly from 1.15–1.23 Å on Pd to 1.60–1.68 Å on PdZn surfaces, similar to the 1.63 Å for the monometallic Zn surface; the C–H distance on Cu surfaces also spans from 1.61–1.63 Å, whilst the distance on the CuPd (110) surface is intermediary to Cu and Pd (1.43 Å). The results imply that the H transfer from the metal surface to the C is most energy intensive for the Pd surface, while the CO₂ activation and bending requires the most energy on Cu, Zn and PdZn surfaces. Further, both the HCOO^b* and HCOO^m* species have stronger (more exothermic) E_ads on lower coordination (and higher energy) surfaces; the HCOO^m* species is also noted as less stable than CO₂^δ− on Pd and CuPd surfaces. The results suggest that the mechanism of CO₂ hydrogenation towards the formate species may be surface dependent: the more exposed H on the PdZn surfaces, relative to Pd, means CO₂ can react with H* directly from the gas phase, i.e., an Eley–Rideal-type mechanism, rather than the Langmuir–Hinshelwood mechanism.

Forming HCOO from physisorbed CO₂ on Pd surfaces has a higher E_a than its decomposition, with E_a(CO₂ + H), 0.33, 0.11 and 0.04 eV greater on Pd (111), (100) and (110), respectively. Cu surfaces show a lower E_a(CO₂ + H), and hence formate remains stable once the carbon is hydrogenated. On Zn (0001), a low E_a(CO₂ + H) of 0.54 eV is calculated assuming a chemisorbed H, but the TS on Zn (0001) is noted as being the highest in energy across all modelled surfaces, with respect to gas substrates, due to the unfavourable Zn–H interaction, visible on the reaction energy profile in Fig. 6. Consideration of the chemisorbed CO₂^δ− species for PdZn(101) gives a low E_a(CO₂^δ− + H) of 0.30 eV, and the barrierless E_a(CO₂^δ− + H) on PdZn (110), indicates that transfer of the H to form a C–H bond requires less energy than on Pd, and the relative E_ads of CO₂ and CO₂^δ− constitutes most of the barrier. We can therefore conclude that the CO₂ activation remains challenging on the surface of the pristine PdZn surfaces, but hydrogenating CO₂ is more efficient than on Pd.

To aid analysis, Pearson correlation was calculated for geometric and energetic quantities related to reactants, products, and TS structures, as shown in Table S4 the SI. The geometries and energies of TS structures are observed to be uncorrelated with the properties of the initial structures in the reactions; however, a positive linear correlation (0.94) of E_a(CO₂ + H) and E_ads of HCOO^m* is observed across Pd (111), (100), (110), Cu (100), (110), CuPd (110), PdZn (101) and (110) surfaces (Fig. S3 in the SI). Simple relationships between adsorbed intermediates and reaction energetics have been previously used successfully to understand and predict catalytic activity for electrochemical CO₂ reduction.^68,69 The trends presented are promising for investigating other Pd-based alloys, where modelling of HCOO^m* could be used to efficiently estimate the E_a(CO₂ + H) without running expensive NEB calculations for transition state searches.

Here, the correlation can be considered meaningful as the HCOO^m* structure captures the interaction of the surface with hydrogen and with the oxygen atoms in the intermediates. The Hammond postulate, coupled with the observed reaction energetics, reasserts that relationship of the HCOO^m* to the TS is sensible, as endergonic reactions have TS structures more akin to the product than the reactants. The surface on which no metastable HCOO^m* was found, Zn (0001), has shown no metastable CO₂^δ− and overall large E_a(CO₂ + H); thus, the absence of the metastable geometry is also informative, as it confirms that the product is unlikely to form. From the trends observed, E_ads(HCOO^m)* emerges as an important target property for further SGD investigation, supporting consideration of catalytically active Pd-based alloys for CO₂ hydrogenation to formate. The chemical context is critical for SGD-supported simulation design and meaningful interpretation of the SGD-derived rules.

3.4 Pd single-atom alloys

To broaden the design options for metal catalysts, the observations in Sections 3.1–3.3 were used as a platform for an SGD analysis of the HCOO^m* adsorbed on SAA Pd surfaces. The investigation focused on Pd–M alloy candidates that facilitate the desired formate pathway of CO₂ hydrogenation, with the correlation between the E_ads(HCOO^m*) and E_a(CO₂ + H) (Fig. S3) allowing us to use E_ads(HCOO^m*) as a target property.

A dataset was prepared by performing a survey of HCOO^m* adsorption across a range of single-atom alloy Pd (111), (100), (110) and (211) surfaces. The results were filtered to ensure that the adsorbate does not scission and the hydrogen does not dissociate from the surface, ensuring that one H–M and at least one O–M bond exists, where M is one of the following species considered: Pd (i.e., non-alloy surface), Co, Cu, Ga, Ir, Ni, Os, Pd, Pt, Rh, Ru, or Zn. The resulting dataset contains 49 symmetrically inequivalent structures, with the lowest E_ads for each SAA shown in Table S5 in the SI.

From the 49 structures considered, 32 result in E_ads(HCOO^m*) > 0 eV (Fig. 7a). To identify subgroups (SGs) of the SAAs that strongly bond HCOO^m*, the normalized negative mean shift utility function was applied when screening data, as described in Section S5 in the SI. Thus the resulting SGD search favours subgroups (SGs) with low mean target values (i.e., exothermic). The SG with 10 data points and desirable average E_ads(HCOO^m*) of −0.30 eV is identified (Fig. 7a, blue bins). The SG includes bridge sites of the (111) and (100) surface, the short and long bridges in the (110) termination, and the bridges located at the step of the (211) surface. The SG includes 3 SAs that are associated with the favourable performance: Co, Ru, and Os. Os is used for homogeneous CO₂ hydrogenation^70,71 but the performance of supported Os heterogeneous nanoparticle catalysts for CO and CO₂ reduction is not strong.^72,73 Ru is a well-known hydrogenation catalyst, and supported Ru nanoparticle catalysts have excelled for CO₂ hydrogenation to formate.⁷⁴ A significant downside of Os and Ru is they are inherently scarce, similar to Pd, and therefore identifying more earth-abundant metals catalysts for CO₂ hydrogenation with the subgroup design rules is required.

Fig. 7 (a) Adsorption energy distribution of the HCOO^m* data set (grey bins) and of the subgroup identified by SGD (blue bins). The bins containing SAs tested after screening for candidate SAAs are shown in orange (Nb, Mo, W) and purple (Ge). (b) Adsorption energy of HCOO^m* on SAAs with respect to the average electron affinity of the metal atoms interacting with the adsorbate (EA_site), one of the key parameters of the selected SG. The highlighted region defines the thresholds given by the SG rules.

The SGD identifies a conjunction (“∧”) of two inequalities that define the selected SG:

EA_SA ≤ 1.192 eV ∧ 0.583 ≤ EA_site ≤ 0.838 eV,

where EA is the electron affinity and the subscript denotes the property belongs to the dopant (SA) or the surface (site). From a chemical perspective EA_SA and EA_site are sensible as they measure the active site’s tendency to withdraw electrons, impacting the charge distribution and its redox flexibility providing means of intermediate stabilisation. Applying the SGD rules allows the search of SAAs suitable for providing strong surface-HCOO^m* binding, which is our focus as it facilitates the Eley–Rideal mechanism and removes the pre-requisite of CO₂ chemisorption for the formate CO₂ hydrogenation pathway.

To further our investigation, the bridge site of Pd (100) was combined with 22 new SA elements for candidate screening. The dopants considered are Sc, Ti, V, Cr, Mn, Fe, Ge, Y, Zr, Nb, Mo, Ag, Cd, In, Sn, Hf, Ta, W, Re, Au, Hg, and Tl. Applying the SG rules reduces the viable SAAs to 4 dopants in this new sample: Cr, Nb, Mo, and W. The validity of the predictions are confirmed by DFT calculations (Table 2), and the E_ads(HCOO^m*) for the SAA with Nb, Mo, and W is very favourable at −0.86, −0.75, and −0.75 eV, respectively (Fig. 7a orange bins). The calculated adsorption energies are below the minimum value in the training data set. To prevent confirmation bias, we also tested two poorly performing dopants not selected by the SGD rules, Ge and Ag. DFT calculations resulted in unfavourable SAA candidates, with the Ge SAA resulting in a highly endothermic metastable structure, with E_ads(HCOO^m*) of 0.70 eV (Fig. 7a purple bin), and no metastable HCOO^m* was observed on the Ag SAA. The DFT calculations confirm the value of the SG rules for identifying dopants capable of enhancing the interaction of HCOO^m* with Pd-based SAAs.

Table 2 Results for the screening of dopants in SAAs with strong binding to HCOO^m*. The rules identified by SGD correctly predict strong interactions for the selected SAAs with Nb, Mo, and W. Moreover, those SAs not selected by the rules, Ge and Ag, do not provide stable HCOO^m* structures. Dashes denote where no minima HCOO^m* structures found after optimisation

Surface	Adsorption site	SA	Selected by SG rules?	Eads(HCOO^m*)/eV
Pd (100)	Bridge	Nb	Yes	−0.86
		Mo		−0.75
		W		−0.75
		Ge	No	0.70
		Ag		—

---

Though Pd-based alloy catalysts with Mo, No and W are not clearly visible in the literature, several studies have explored Pd catalysts on Mo-, Nb- and W-based materials for CO₂ hydrogenation in various types of catalytic CO₂ hydrogenation process. In terms of thermocatalytic CO₂ hydrogenation, Mo₂C-supported Pd nanoparticles show significantly enhanced activity, with lower apparent activation energies for CO₂ hydrogenation to formate than their Ru–Pd/Mo₂C counterparts (where Ru merely stabilises the Pd NPs).⁷⁵ Likewise, Pd/MoO₃ catalysts show substantially higher catalytic activity than Cu/ZnO catalysts in liquid-phase CO₂ hydrogenation under mild conditions.⁷⁶ Pd/NbC and Pd/WC materials have been reported as potent CO₂ electroreduction catalysts, owing to favourable Pd–M electronic interactions.⁷⁷ In thermo-photo-catalytic reverse water–gas shift reaction (RWGS), small Pd NPs on Nb₂O₅ are highly selective towards CO, suppressing CH₄ formation typical of Pd catalysts.⁷⁸ Finally, Nb-doped Pd/TiO₂ also demonstrates synergistic thermo-photo-catalytic production of syngas (CO + H₂) from methanol/water mixtures, with notable effects even at 0.5 wt% loadings.⁷⁹

In summary, there are examples where the Mo, Nb and W atoms may play an important role for CO₂ hydrogenation; Mo–Pd, Nb–Pd, W–Pd alloy nanoparticle-based catalysts for CO₂ reduction are not well explored, and the metal oxide- or metal carbide-supported Pd catalysts reported thus far are not typically thoroughly tested for the presence of these alloys. Only the PdMo alloy catalysts have been reported, with ability to hydrogenate CO₂ to methanol at 25 °C and 9 bar pressure, which are much milder conditions than the conventional CZA catalyst, though at lower turnover rates.⁸⁰ Computational works show that the MoNb-doped Pd surfaces with dual single-atom sites are predicted to be especially active for CO₂ activation⁴² and near-surface Pd/W alloys favour the formation of HCOO* over COOH* in electrocatalytic CO₂ reduction.⁸¹

In summary, the subgroup discovery has been successfully applied to identify promising Pd–M alloy candidates for CO₂ hydrogenation to methanol. Mo, Nb and W are highlighted as potent alloy candidates enhancing the selectivity of Pd catalysts by facilitating the Eley–Rideal mechanism of CO₂ hydrogenation to formate. The mechanistic change alleviates the challenge of CO₂ pre-activation and allows the reaction to proceed via the formate pathway of methanol synthesis. The Pd–M alloy SGD analysis presented herein provides strong evidence of enhanced CO₂ activation towards formate but is not informed on further mechanistic aspects of the reaction, which could ultimately dictate the selectivity of CO₂ hydrogenation on these alloy catalysts. The present calculations neglect several catalytically relevant factors, including the influence of temperature and pressure on surface composition and reaction pathways; future work should therefore explicitly account for these effects. Future first-principles studies are required to confirm the viability of Mo–Pd, Nb–Pd, Pd–W alloy formation and stability under reaction conditions, and to explore further the mechanistic aspects of CO₂ hydrogenation on these materials.

4 Conclusions

The synthesis of methanol from CO₂ and hydrogen using Pd-based alloy catalysts holds great potential for sustainable energy storage. Gaining deeper knowledge of the initial stages of CO₂ hydrogenation on Pd-based catalysts is essential for understanding the variation in reactivity observed upon alloying, which is important for the strategic design of novel catalysts. Here, we have investigated the initial CO₂ hydrogenation to formate across monometallic FCC Pd (111), (100), (110), FCC Cu (111), (100), (110), HCP Zn (0001) surfaces and BCC CuPd (110) and BCT PdZn (101) and (110) alloy facets. Cu and PdZn surfaces have the lowest hydrogenation barriers. CuPd (110) has a similar activation energy to Pd surfaces, but the formate product is much less stable on the surface. CO₂ activation constitutes a major part of the CO₂ hydrogenation barrier on PdZn and Cu surfaces, thus making the chemical pre-adsorption of CO₂ unlikely and suggests an Eley–Rideal mechanism as opposed to the Langmuir–Hinshelwood model. Statistical analysis shows that there are limited correlations between reactants and transition state structures, but the adsorption energy of the HCOO^m product does correlate with the activation energy required for hydrogenation of the CO₂ species.

Subgroup discovery was applied on a dataset derived from a HCOO^m* adsorption survey on a range of SAA Pd (111), (100), (110) and (211) surfaces, doped with Co, Cu, Ga, Ir, Ni, Os, Pd, Pt, Rh, Ru or Zn. The chosen 13 candidate features relating to the dopant or the adsorption site of HCOO^m* were used to capture the SAA surfaces’ electronic and geometric properties. The combination of DFT and SGD analysis shows that the electron affinity (EA) of the dopant, which is an intrinsic material property, can be used to predict which SAA surfaces facilitate HCOO^m* adsorption. The resulting selectivity criteria allows us to identify Nb, Mo, and W as candidates to alloy with Pd catalysts and significantly improve the capability of CO₂ hydrogenation to formate. The materials have the potential to enable the formate pathway of methanol synthesis on Pd-based catalysts with much improved precious metal utilisation.

The workflow and results from this investigation can lead the way for theory-led rational design of alloy catalysts, in this case specifically targeting mechanistic aspects of CO₂ hydrogenation. The general workflow offers potential applicability to other reactions making use of alloy catalysts, accelerating catalyst discovery.

Author contributions

Igor Kowalec: conceptualisation, data curation, formal analysis, investigation, methodology, resources, validation, visualisation, roles/writing – original draft; Herzain Rivera Arrieta: data curation, formal analysis, methodology, resources, software, roles/writing – original draft; Zhongwei Lu: data curation, formal analysis, roles/writing – original draft; Lucas Foppa: formal analysis, methodology, software, supervision, writing – review & editing; Matthias Scheffler: formal analysis, funding acquisition, methodology, software, supervision, writing – review & editing; Richard Catlow: conceptualisation, funding acquisition, methodology, supervision, writing – review & editing; Andrew Logsdail: conceptualisation, formal analysis, funding acquisition, methodology, project administration, resources, supervision, writing – review & editing.

Conflicts of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

The inputs and output files of all electronic structure calculations have been uploaded as a dataset to the NOMAD repository at DOI: https://doi.org/10.17172/NOMAD/2026.02.02-1.

Supplementary information (SI) is available. See DOI: https://doi.org/10.1039/d5fd00125k.

Acknowledgements

The authors are grateful for funding by the EPSRC Centre-to Centre Project (Grant reference: EP/S030468/1). We are grateful to Graham Hutchings and David Willock for useful discussions. AJL acknowledges funding by the UKRI Future Leaders Fellowship program (MR/T018372/1, MR/Y034279/1). ZL acknowledges Cardiff University-China Scholarship Council Research Excellence Scholarship. The authors acknowledge computational resources and support from: the Supercomputing Wales project, which is part-funded by the European Regional Development Fund (ERDF) via the Welsh Government; the UK National Supercomputing Service YOUNG, accessed via membership of the Materials Chemistry Consortium which is funded by Engineering and Physical Sciences Research Council (EP/L000202, EP/R029431, EP/T022213, EP/X035859); and the Isambard UK National Tier-2 HPC Service operated by GW4 and the UK Met Office, and funded by EPSRC (EP/P020224/1, EP/T022078/1). This work was partially funded by the ERC Advanced Grant TEC1p (European Research Council, Grant Agreement No. 740233). The research was performed within the framework of the Max Planck-Cardiff Centre on the Fundamentals of Heterogeneous Catalysis (FUNCAT) at Cardiff University.

References

1. D. A. Jurado A., M. D. Higham, Y. R. Poh, C. R. A. Catlow and I. Krossing, J. Catal., 2025, 446, 115997.
2. A. Prašnikar, A. Pavlišič, F. Ruiz-Zepeda, J. Kovač and B. Likozar, Ind. Eng. Chem. Res., 2019, 58, 13021–13029.
3. W. Wang, S. Wang, X. Ma and J. Gong, Chem. Soc. Rev., 2011, 40, 3703–3727.
4. M. D. Porosoff, B. Yan and J. G. Chen, Energy Environ. Sci., 2016, 9, 62–73.
5. A. L. Maulana, R. I. D. Putra, A. G. Saputro, M. K. Agusta, Nugraha and H. K. Dipojono, Phys. Chem. Chem. Phys., 2019, 21, 20276–20286.
6. M. Zhang, Y. Wu, M. Dou and Y. Yu, Catal. Lett., 2018, 148, 2935–2944.
7. M. D. Higham, M. G. Quesne and C. R. A. Catlow, Dalton Trans., 2020, 49, 8478–8497.
8. I. Kowalec, L. Kabalan, C. R. A. Catlow and A. J. Logsdail, Phys. Chem. Chem. Phys., 2022, 24, 9360–9373.
9. F. Brix, V. Desbuis, L. Piccolo and É. Gaudry, J. Phys. Chem. Lett., 2020, 11, 7672–7678.
10. L. Kabalan, I. Kowalec, C. R. A. Catlow and A. J. Logsdail, Phys. Chem. Chem. Phys., 2021, 23, 14649–14661.
11. K. P. Kuhl, T. Hatsukade, E. R. Cave, D. N. Abram, J. Kibsgaard and T. F. Jaramillo, J. Am. Chem. Soc., 2014, 136, 14107–14113.
12. M. W. Roberts, J. Mol. Catal., 1992, 74, 11–22.
13. P. R. Davies and M. W. Roberts, J. Chem. Soc. Faraday Trans., 1992, 88, 361–368.
14. A. F. Carley, S. Yan and M. W. Roberts, J. Chem. Soc. Faraday Trans., 1990, 86, 2701–2710.
15. K. Blaziak, D. Tzeli, S. S. Xantheas and E. Uggerud, Phys. Chem. Chem. Phys., 2018, 20, 25495–25505.
16. T.-W. Jiang, X. Qin, K. Ye, W.-Y. Zhang, H. Li, W. Liu, S. Huo, X.-G. Zhang, K. Jiang and W.-B. Cai, Appl. Catal., B, 2023, 334, 122815.
17. L. C. Grabow and M. Mavrikakis, ACS Catal., 2011, 1, 365–384.
18. M. Huš, D. Kopač, N. S. Štefančič, D. L. Jurković, V. D. B. C. Dasireddy and B. Likozar, Catal. Sci. Technol., 2017, 7, 5900–5913.
19. A. J. Medford, J. Sehested, J. Rossmeisl, I. Chorkendorff, F. Studt, J. K. Nørskov and P. G. Moses, J. Catal., 2014, 309, 397–407.
20. J. Díez-Ramírez, P. Sánchez, J. L. Valverde and F. Dorado, J. CO₂ Util., 2016, 16, 375–383.
21. J. Díez-Ramírez, J. L. Valverde, P. Sánchez and F. Dorado, Catal. Lett., 2016, 146, 373–382.
22. E. J. Peterson, B. Halevi, B. Kiefer, M. N. Spilde, A. K. Datye, J. Peterson, L. Daemen, A. Llobet and H. Nakotte, J. Alloys Compd., 2011, 509, 1463–1470.
23. M. Miyazaki, S. Furukawa, T. Takayama, S. Yamazoe and T. Komatsu, ACS Appl. Nano Mater., 2019, 2, 3307–3314.
24. M. Heggen, S. Penner, M. Friedrich, R. E. Dunin-Borkowski and M. Armbrüster, J. Phys. Chem. C, 2016, 120, 10460–10465.
25. M. Friedrich, S. Penner, M. Heggen and M. Armbrüster, Angew. Chem., Int. Ed., 2013, 52, 4389–4392.
26. E. Nowicka, S. M. Althahban, Y. Luo, R. Kriegel, G. Shaw, D. J. Morgan, Q. He, M. Watanabe, M. Armbrüster, C. J. Kiely and G. J. Hutchings, Catal. Sci. Technol., 2018, 8, 5848–5857.
27. R. Qiu, Z. Ding, Y. Xu, Q. Yang, K. Sun and R. Hou, Appl. Surf. Sci., 2021, 544, 148974.
28. L. Kabalan, I. Kowalec, S. Rigamonti, M. Troppenz, C. Draxl, C. R. A. Catlow and A. J. Logsdail, J. Phys.: Condens. Matter, 2023, 35, 405402.
29. S. Ma, M. Sadakiyo, M. Heima, R. Luo, R. T. Haasch, J. I. Gold, M. Yamauchi and P. J. A. Kenis, J. Am. Chem. Soc., 2017, 139, 47–50.
30. L. Zhu, Y. Lin, K. Liu, E. Cortés, H. Li, J. Hu, A. Yamaguchi, X. Liu, M. Miyauchi, J. Fu and M. Liu, Chin. J. Catal., 2021, 42, 1500–1508.
31. I. Melián-Cabrera, M. L. Granados and J. L. G. Fierro, J. Catal., 2002, 210, 285–294.
32. V. S. S. Mosali, L. Li, G. Puxty, M. D. Horne, A. M. Bond and J. Zhang, ChemElectroChem, 2022, 9, e202101227.
33. W. J. Wang, S. Hwang, T. Kim, S. Ha and L. Scudiero, Electrochim. Acta, 2021, 387, 138531.
34. H. Pan, B. Ma, L. Zhou, Y. Hu, M. Shakouri, Y. Guo, X. Liu and Y. Wang, ACS Sustain. Chem. Eng., 2023, 11, 7489–7499.
35. M. Bowker, ChemCatChem, 2019, 11, 4238–4246.
36. M. B. Boucher, B. Zugic, G. Cladaras, J. Kammert, M. D. Marcinkowski, T. J. Lawton, E. C. H. Sykes and M. Flytzani-Stephanopoulos, Phys. Chem. Chem. Phys., 2013, 15, 12187–12196.
37. G. Giannakakis, M. Flytzani-Stephanopoulos and E. C. H. Sykes, Acc. Chem. Res., 2019, 52, 237–247.
38. R. T. Hannagan, G. Giannakakis, M. Flytzani-Stephanopoulos and E. C. H. Sykes, Chem. Rev., 2020, 120, 12044–12088.
39. J. Schumann, M. Stamatakis, A. Michaelides and R. Réocreux, Nat. Chem., 2024, 16, 749–754.
40. Y. Wang, M. Zheng, X. Zhou, Q. Pan and M. Li, Molecules, 2023, 28, 4606.
41. F. J. Caparrós, L. Soler, M. D. Rossell, I. Angurell, L. Piccolo, O. Rossell and J. Llorca, ChemCatChem, 2018, 10, 2365–2369.
42. H. I. Rivera-Arrieta and L. Foppa, ACS Catal., 2025, 15, 2916–2926.
43. S. Wrobel, in Principles of Data Mining and Knowledge Discovery, ed. J. Komorowski and J. Zytkow, Springer, Berlin, Heidelberg, 1997, pp. 78–87.
44. V. Blum, R. Gehrke, F. Hanke, P. Havu, V. Havu, X. Ren, K. Reuter and M. Scheffler, Comput. Phys. Commun., 2009, 180, 2175–2196.
45. A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dulak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. B. Jensen, J. Kermode, J. R. Kitchin, E. L. Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. B. Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiøtz, O. Schütt, M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. Zeng and K. W. Jacobsen, J. Phys.: Condens. Matter, 2017, 29, 273002.
46. E. van Lenthe, J. G. Snijders and E. J. Baerends, J. Chem. Phys., 1996, 105, 6505–6516.
47. J. Wellendorff, K. T. Lundgaard, K. W. Jacobsen and T. Bligaard, J. Chem. Phys., 2014, 140, 144107.
48. S. Lehtola, C. Steigemann, M. J. T. Oliveira and M. A. L. Marques, SoftwareX, 2018, 7, 1–5.
49. J. S. Gómez-Jeria, J. Chil. Chem. Soc., 2009, 54, 482–485.
50. M. Yamauchi and T. Tsukuda, Dalton Trans., 2011, 40, 4842–4845.
51. J. W. M. Crawley, I. E. Gow, N. Lawes, I. Kowalec, L. Kabalan, C. R. A. Catlow, A. J. Logsdail, S. H. Taylor, N. F. Dummer and G. J. Hutchings, Chem. Rev., 2022, 122, 6795–6849.
52. M. H. Hansen, J. A. G. Torres, P. C. Jennings, Z. Wang, J. R. Boes, O. G. Mamun and T. Bligaard, arXiv, 2019, preprint, arXiv:1904.00904,  DOI:10.48550/arXiv.1904.00904.
53. G. Henkelman, B. P. Uberuaga and H. Jónsson, J. Chem. Phys., 2000, 113, 9901–9904.
54. H. I. Rivera-Arrieta, NOMAD Repository, 2024,  DOI:10.17172/NOMAD/2024.05.10-1.
55. I. Batatia, P. Benner, Y. Chiang, A. M. Elena, D. P. Kovács, J. Riebesell, X. R. Advincula, M. Asta, M. Avaylon, W. J. Baldwin, F. Berger, N. Bernstein, A. Bhowmik, S. M. Blau, V. Cărare, J. P. Darby, S. De, F. Della Pia, V. L. Deringer, R. Elijošius, Z. El-Machachi, F. Falcioni, E. Fako, A. C. Ferrari, A. Genreith-Schriever, J. George, R. E. A. Goodall, C. P. Grey, P. Grigorev, S. Han, W. Handley, H. H. Heenen, K. Hermansson, C. Holm, J. Jaafar, S. Hofmann, K. S. Jakob, H. Jung, V. Kapil, A. D. Kaplan, N. Karimitari, J. R. Kermode, N. Kroupa, J. Kullgren, M. C. Kuner, D. Kuryla, G. Liepuoniute, J. T. Margraf, I.-B. Magdău, A. Michaelides, J. H. Moore, A. A. Naik, S. P. Niblett, S. W. Norwood, N. O’Neill, C. Ortner, K. A. Persson, K. Reuter, A. S. Rosen, L. L. Schaaf, C. Schran, B. X. Shi, E. Sivonxay, T. K. Stenczel, V. Svahn, C. Sutton, T. D. Swinburne, J. Tilly, C. van der Oord, E. Varga-Umbrich, T. Vegge, M. Vondrák, Y. Wang, W. C. Witt, F. Zills and G. Csányi, arXiv, 2024, preprint, arXiv:2401.00096,  DOI:10.48550/arXiv.2401.00096.
56. S. Takamoto, C. Shinagawa, D. Motoki, K. Nakago, W. Li, I. Kurata, T. Watanabe, Y. Yayama, H. Iriguchi, Y. Asano, T. Onodera, T. Ishii, T. Kudo, H. Ono, R. Sawada, R. Ishitani, M. Ong, T. Yamaguchi, T. Kataoka, A. Hayashi, N. Charoenphakdee and T. Ibuka, Nat. Commun., 2022, 13, 2991.
57. S. Grimme, J. Comput. Chem., 2004, 25, 1463–1473.
58. B. R. Goldsmith, M. Boley, J. Vreeken, M. Scheffler and L. M. Ghiringhelli, New J. Phys., 2017, 19, 013031.
59. L. Foppa, C. Sutton, L. M. Ghiringhelli, S. De, P. Löser, S. A. Schunk, A. Schäfer and M. Scheffler, ACS Catal., 2022, 12, 2223–2232.
60. P. S. Schmidt and K. S. Thygesen, J. Phys. Chem. C, 2018, 122, 4381–4390.
61. J. Quan, F. Muttaqien, T. Kondo, T. Kozarashi, T. Mogi, T. Imabayashi, Y. Hamamoto, K. Inagaki, I. Hamada, Y. Morikawa and J. Nakamura, Nat. Chem., 2019, 11, 722–729.
62. E. Iglesia, J. Phys. Chem., 1986, 90, 5272–5274.
63. G. J. Millar, C. H. Rochester and K. C. Waugh, J. Chem. Soc. Faraday Trans., 1991, 87, 1491–1496.
64. Y. Yang, C. A. Mims, R. S. Disselkamp, J. H. Kwak, C. H. F. Peden and C. T. Campbell, J. Phys. Chem. C, 2010, 114, 17205–17211.
65. J. Quan, T. Kondo, G. Wang and J. Nakamura, Angew. Chem., Int. Ed., 2017, 56, 3496–3500.
66. G. Wang, Y. Morikawa, T. Matsumoto and J. Nakamura, J. Phys. Chem. B, 2006, 110, 9–11.
67. S. Li, S. Rangarajan, J. Scaranto and M. Mavrikakis, Surf. Sci., 2021, 709, 121846.
68. H. A. Hansen, J. B. Varley, A. A. Peterson and J. K. Nørskov, J. Phys. Chem. Lett., 2013, 4, 388–392.
69. X. Liu, J. Xiao, H. Peng, X. Hong, K. Chan and J. K. Nørskov, Nat. Commun., 2017, 8, 15438.
70. S. S. P. Malaza and B. C. E. Makhubela, J. CO₂ Util., 2020, 39, 101149.
71. H. W. Choi and E. L. Muetterties, Inorg. Chem., 1981, 20, 2664–2667.
72. S. D. Jackson, R. B. Moyes, P. B. Wells and R. Whyman, J. Catal., 1984, 86, 342–358.
73. H. Knözinger, Y. Zhao, B. Tesche, R. Barth, R. Epstein, B. C. Gates and J. P. Scott, Faraday Discuss. Chem. Soc., 1981, 72, 53–71.
74. S. J. Louis Anandaraj, L. Kang, S. DeBeer, A. Bordet and W. Leitner, Small, 2023, 19, 2206806.
75. C. E. Mitchell, U. Terranova, I. Alshibane, D. J. Morgan, T. E. Davies, Q. He, J. S. J. Hargreaves, M. Sankar and N. H. de Leeuw, New J. Chem., 2019, 43, 13985–13997.
76. Y. Kuwahara, T. Mihogi, K. Hamahara, K. Kusu, H. Kobayashi and H. Yamashita, Chem. Sci., 2021, 12, 9902–9915.
77. J. Wang, S. Kattel, C. J. Hawxhurst, J. H. Lee, B. M. Tackett, K. Chang, N. Rui, C.-J. Liu and J. G. Chen, Angew. Chem., Int. Ed., 2019, 58, 6271–6275.
78. J. Jia, P. G. O’Brien, L. He, Q. Qiao, T. Fei, L. M. Reyes, T. E. Burrow, Y. Dong, K. Liao, M. Varela, S. J. Pennycook, M. Hmadeh, A. S. Helmy, N. P. Kherani, D. D. Perovic and G. A. Ozin, Adv. Sci., 2016, 3, 1600189.
79. U. Caudillo-Flores, R. Sayago, A. Ares-Dorado, S. Fuentes-Moyado, M. Fernández-García and A. Kubacka, ACS Sustain. Chem. Eng., 2023, 11, 3896–3906.
80. H. Sugiyama, M. Miyazaki, M. Sasase, M. Kitano and H. Hosono, J. Am. Chem. Soc., 2023, 145, 9410–9416.
81. Z. Zhao and G. Lu, ACS Catal., 2018, 8, 3885–3894.

---