Catalytic Potential of Post-Transition Metal Doped Graphene-Based Single-Atom Catalysts for the CO2 Electroreduction Reaction.

Catalysts are required to ensure electrochemical reduction of CO2 to fuels proceeds at industrially acceptable rates and yields. As such, highly active and selective catalysts must be developed. Herein, a density functional theory study of p-block element and noble metal doped graphene-based single-atom catalysts in two defect sites for the electrochemical reduction of CO2 to CO and HCOOH is systematically undertaken. It is found that on all of the systems considered, the thermodynamic product is HCOOH. Pb/C3 , Pb/N4 and Sn/C3 are identified as having the lowest overpotential for HCOOH production while Al/C3 , Al/N4 , Au/C3 and Ga/C3 are identified as having the potential to form higher order products due to the strength of binding of adsorbed HCOOH.

Introduction

Global climate change is driven by the over‐abundance of CO2 in the atmosphere. One mitigation strategy is to convert atmospheric CO2 to fuel sources via the CO2 electroreduction reaction (CO2RR). However, the C=O bond is particularly stable (∼799 kJ mol−1) and, therefore, a catalyst is required to enable this reaction to proceed at appreciable rates and yields. Furthermore, the CO2RR is a highly complex process which results in a wide array of reaction products.

The mechanism of reaction on some catalysts is initiated by CO2 physisorbing to the surface, before a transfer of partial charge to the antibonding orbital of physisorbed CO2 causes it to bend, thus rendering it ‘activated.’[ , , ] In this pathway, the bent chemisorbed CO2 is then reacted with a H to form adsorbed COOH. From here, a multitude of different products can form, including CO, CH4, C2H4 and C2H5OH.[ , , ] However, if the catalyst has a particularly strong affinity for O binding, the reaction can proceed without the CO2 molecule bending.[ , ] In this reaction pathway, the adsorbed CO2 reacts with H to give OCHO, the key intermediate for formate (HCOO−) or formic acid (HCOOH) production. These competing mechanistic pathways make designing an effective catalyst for this reaction a non‐trivial task, because the catalyst must be active for the CO2RR and also highly selective for a desired product.

Many different catalysts have been developed for the CO2RR. The simplest of these are monometallic catalysts. However, Cu is the only monometallic transition metal capable of performing the CO2RR to produce alcohols and hydrocarbons at reasonable (∼70 %) efficiencies.[ , , , ] Bimetallic catalysts and single‐atom alloy catalysts have also been considered, aiming to couple the selectivity of one metal to the activity of another for optimal catalyst performance.[ , ] In addition to mixing metals to improve catalytic ability, different surface features tune the activity and selectivity of these materials.[ , ] It is well known that planar sites are comparably unreactive, while considerably higher activity can be accessed on step, kink, edge and terrace sites.[ , , ]

In recent years, graphene‐based systems have been functionalised with metal and non‐metal dopants to create so‐called ‘single‐atom’ catalysts (SACs). SACs are unique because, at the atomic size, dopant atoms can have properties quite distinct from the parent material at the nanometer or bulk size regimes thus providing a new avenue of investigation for catalyst candidates. It is intuitive that different dopant species in SACs would exhibit different functionality for the CO2RR, as confirmed by numerous studies.[ , , , , ] However, changing the metal dopant is not the only way to alter the catalytic efficiency of these materials. Co‐doping of metal and non‐metal species into defect sites can result in synergistic effects that alter the behaviour of these materials further. It has been shown that decreasing the number of coordinating N atoms can increase a system's Faradaic efficiency for the CO2RR.[ , , , , , ] In SACs, the dominant product formed is CO due to the single‐atom sites being sufficiently far apart that C−C bonds are unable to form. However, it has been shown that the proximity of sites in Cu SACs allows fine‐tuning of the material to produce a variety of hydrocarbons. The tuning of metal dopant, defect site structure, co‐doping, proximity of active sites and changing the support material provides many combinations of materials available for optimisation.

Within these classes of catalysts (monometallic, bimetallic and SACs) many different arrangements and combinations of metals and non‐metals have been trailed in an effort to develop the optimal catalyst.[ , , , , ] However, this is effectively a trial‐and‐error method to catalyst design and is an extremely laborious process. To circumvent this, relationships have been recognised to predict catalytic activity and, thus, identify particularly promising catalytic materials.

The best known guiding principles in the prediction of catalytic promise are scaling relations and Brønsted‐Evans‐Polyani (BEP) relations. Scaling relations state that the adsorption of a simple molecule correlates to the adsorption energy of key intermediates in a reaction pathway when bound to the surface through the same element. The manipulation of scaling relations can result in the formation of so‐called limiting potential volcano plots which are able to not only predict the rate‐limiting step of a reaction pathway, but also the best candidate material for a particular reaction. Scaling relations have been used extensively on mono‐ and bimetallic catalysts.[ , , ] However, they have been met with mixed success for SACs with some studies finding scaling relations to be adhered to while other studies have found that SACs break scaling relations.[ , ]

BEP relations allow the transition state energy to be inferred from the adsorption energy of intermediates. While the energetic barrier for a reaction can only be accurately determined by locating the transition state, BEP relationships suggest that the reaction energy is proportional to the activation energy barrier.[ , ] If the reaction energy is exergonic, based on the BEP relationships, this will likely lead to a low transition state energy. BEP relations have been shown to hold for the CO2RR on bimetallic alloys, Pb surfaces, and dilute single‐atom alloys but the viability of BEP relations are yet to be examined explicitly for graphene‐based SACs.

To aid in identifying materials with catalytic potential, activity descriptors have been developed to determine the binding energy of key intermediates. In conjunction with scaling and BEP relations, activity descriptors can be used to describe and predict the catalytic activity of different materials. The d‐band centre is the most widely known descriptor with the average position of valence d‐shell with respect to the Fermi level being correlated to the binding energy of reaction intermediates for transition metals.[ , ] The d‐band centre model has been widely used as a descriptor for monometallic systems,[ , , ] however, using this parameter for graphene‐based transition‐metal doped SACs is somewhat complicated because the metal d‐orbitals interact with the C p‐orbitals of the graphene.

Other activity descriptors have also been proposed that pertain specifically to SACs, including Bader charge, vacancy formation energy, difference in binding energy of adsorbed CO and adsorbed H and group number. Perhaps most interestingly, however, is the development of so‐called ‘universal descriptors’,[ , , ] where ‘universal’ implies that the descriptor can account for different coordination environments and metals while correctly predicting activity and selectivity of potential catalytic materials using this single parameter. Descriptors used across a wide variety of catalytic materials provide a rationale for the differing behaviour of catalytic systems and suggest a framework for optimising catalysts.

While descriptors are highly effective, most are centred around the partial filling of the d‐band; the defining feature of the transition‐metals, for which they have been designed to describe. The partially filled d‐band is involved in binding adsorbates and, therefore, is one of the primary influencing factors in determining a transition metal's activity and/or selectivity. However, due to the reliance of these descriptors on the partially filled d‐band they are not readily applicable for non‐transition‐metal doped systems. Nonetheless, p‐block elements appear in the literature as promising catalytic materials, in various forms. For example, an Al metal organic framework was recently found to favour the CO2RR and suppress the hydrogen evolution reaction (HER). Sn, In, Bi and Pb surfaces are selective for HCOOH production with a Faradaic efficiency of ∼90 %.[ , , , , ] In, Sn and Bi can be used to build SACs, with HCOO‐ as their most commonly produced product,[ , , , ] but with specific coordination environments are able to be tailored to produce CO as the major product.[ , ] Thus, it is clear that p‐block elements should be meticulously explored for their catalytic potential. Despite this, screening studies of the p‐block elements for catalysis are non‐existent.

Here, a systematic study is carried out on a range of p‐block element (Al, Bi, Ga, In, Sn, Pb) and Au‐doped graphene‐based SACs in the C3‐site and N4‐site for their potential for the CO2RR. The intermediates for the CO2RR two‐electron transfer products; CO and HCOOH, are examined in detail. Furthermore, a range of possible descriptors for these materials are explored.

Computational Details

Calculation Details

Calculations were carried out in the Vienna ab initio Simulation Package (VASP) using plane wave density functional theory (DFT) with the PBE exchange‐correlational functional and projector‐augmented wave (PAW) method. An energy cutoff of 600 eV was used. Van der Waals interactions were included through the D3 method with Becke‐Johnson damping.[ , ] Calculations were electronically converged to  eV, while ionic convergence was set to 1×10−2 eV to ensure minor optimisation fluctuations in the graphene sheet were avoided. Monkhorst‐Pack k‐point sampling of was selected. Higher k‐points were tested and energies altered by less than 2×10−4 eV atom−1 as the k‐point grid increased (Figure S1).

PAW datasets that included d‐electrons in the core were used for Ga, Bi and Sn. For Au, In and Pb, d‐electrons were included in the valence. For Au, this was due to the d‐electrons being required for bonding. However, for In and Pb, defining the d‐electrons in the valence was necessary to avoid numerical instabilities in the calculations. Careful testing was undertaken on those dopants that could accommodate both d‐electrons within the core and in the valence in the computational setup; Bi, Ga and Sn. The energetics of an adsorbed CO2 were altered by less than 0.02 eV (Table S2) and the Bader charge of the doped site was effectively unchanged upon inclusion of the d‐electrons in the valence (Table S3). Thus, we are confident that all systems are able to be directly compared to one another regardless of the valency of the PAW dataset for the p‐block elements.

Model Setup

The objective of this work was to examine the potential of selected metal dopants in two defect sites in graphene‐based SACs for the CO2RR. The dopants considered in this study were Al, Au, Bi, Ga, In, Pb and Sn. The two sites considered were a graphitic site, where a single C atom was removed from the graphene sheet (hereafter referred to as the C3‐site; Figure S1) and a pyridinic site where a double C vacancy was created and then reacted with N to form a pyridinic‐type arrangement (hereafter referred to as the N4‐site; Figure S2). Unit cell sizes of 4×4 hexagonal C rings were used, resulting in 50 atoms in the C3‐site system and 49 atoms in the N4‐site system, including the metal dopant. Repeating units in the z‐dimension were separated by a minimum of 15 Å of vacuum. Example position files for Al/C3 and Al/N4 can be found in Supplementary Note 2. The stability of the majority of the single‐atom catalysts is examined computationally in a previous work, however, in an experimental context, co‐doping the metal‐dopant with other elements may also help with stabilisation and durability of the active site.[ , ]

In addition to CO2, intermediates for the two‐electron transfer products of electrochemical CO2RR (CO and HCOOH) were adsorbed on the surface, as shown in the balanced half equations 1–4b. An asterisk represents a surface site and an asterisk next to an adsorbate indicates that the adsorbate is adsorbed to the surface site. Initial geometries trialled can be found in Tables S4 and S5.

Method

Here, the computational hydrogen electrode (CHE) was employed, as outlined by Nørskov et al. which accounts for the energy of a proton‐electron pair (H++e−) in an aqueous solution, by assuming:

and that the free energy change of *A+H++e−→*AH is equivalent to *A+ H2(g)→*AH.

As a correction to the electronic energy, Gibbs free energies of all systems are reported here. The Gibbs free energy (G) of reaction intermediates was calculated using the following equation:

where E is the electronic energy, ZPE is the zero point energy, T is the temperature (298.15 K) and S is the entropic contribution from the vibrational motion of the adsorbate. ZPE and S contributions were calculated using harmonic vibrational frequencies obtained from normal mode analysis for adsorbates on the Al/C3‐site and applied to all systems (Table S6). In the calculation of ZPE and TS, only contributions from the adsorbate were considered; contributions from the graphene flake were excluded. Values calculated in the present study are in good agreement with those used in other works.[ , , ] The error introduced as a result of applying the Al/C3 ZPE‐TS correction for each adsorbate on all systems was estimated by carrying out several careful checks (Table S7). The largest deviation in the calculation of ZPE‐TS between the Al/C3 system and all other systems for which the normal mode analysis was carried out was 0.20 eV, meaning that this is a good proxy for maximal error within the current calculations arising from the ZPE and TS corrections that have been applied.

Enthalpic contributions and solvation corrections were not included in this study. Enthalpic contributions have been shown to contribute no more than 0.11 eV and are remarkably consistent between adsorbates (0.08 eV–0.11 eV), resulting in, at most, a negligible 0.03 eV change in relative energetics. Solvation corrections were not included as explicit calculation of a water layer is much too computationally expensive. However, it has previously been shown that the presence of water can stabilise some intermediates; by 0.5 eV for *OH, 0.25 eV for an *R‐OH group (e. g. *COOH) and 0.1 eV for *CO. However solvation corrections are not presented in the literature for *OCHO, thus it is unjustifiable to only apply solvation corrections to selected intermediates.

Free energies were calculated with reference to G(H), G(O) and G(C) and the electronic energy of the optimised clean metal surface (*). These reference states were calculated in keeping with those presented in the literature[ , , ] as:

where n is the number of atoms in the graphene unit cell. For the reference states, ZPE energies were calculated explicitly from the normal mode analysis of vibrational frequency calculations but standard molar entropic contributions, S, were taken from thermodynamic tables to allow for inclusion of vibrational and rotational contributions to the entropy (Table S8). Additionally, free energy diagrams are presented relative to G(*+CO2), which itself is calculated with respect to the reference states as above, to allow for direct comparison of the energy of gaseous and adsorbed species across the free energy diagrams.

To calculate the change in free energy at each step we employ the equation:

where n is the number of electrons transferred in the step, e is the charge on an electron, U is the applied potential and is the free energy in the absence of potential. Here, the applied potential is 0.0 V, therefore is equivalent to . All potentials presented here are vs. the reversible hydrogen electrode (RHE) at pH=0. Full calculations of the free energy of intermediates for the electrochemical CO2RR to the two electron transfer products on the Al/C3 catalyst can be found in Supplementary Note S5.

The theoretical limiting potential (U), specifically at U=0.0 V, is defined as:

thus, the most negative U for each reaction pathway represents the limiting potential for a particular product, on a specific catalyst. This only applies to the electrochemical steps where a are transferred (equations 2a, 2b, 3a, 3b). The adsorption of CO2 (equation 1) and desorption of the product (equations 4a and 4b) are not electrochemical steps, therefore the change in free energy of these steps must be sufficiently small to be able to be overcome thermally under reaction conditions.

Results and Discussion

Competing CO2RR vs HER

The HER is competitive with the CO2RR under electrochemical reduction conditions. Catalysts can be screened by plotting ΔG(*H), the intermediate species in the HER, against possible intermediates in the CO2RR, ΔG(*COOH) and ΔG(*OCHO) (Figure 1). When the free energies of adsorption of CO2RR intermediates are more stabilised than the reaction intermediate of the HER, the catalyst is selective for the CO2RR. If the *H intermediate is more stabilised, the catalyst is likely to be rendered inactive as it is effectively blocked by *H. From Figure 1, those catalysts that are selective for the CO2RR are Al/C3, Al/N4, Pb/C3, Pb/N4, Au/N4, Bi/C3 and Bi/N4. However, here we highlight that we have not included solvation effects in the calculation of the ΔG(*COOH), which would likely increase the competitiveness of doped SACs for the CO2RR, over the HER by stabilising the *COOH intermediate by ∼0.25 eV.[ , ]

Figure 1

Relation between ΔG(*H) and key intermediates in the CO2RR, ΔG(*OCHO) and ΔG(*COOH). Single‐atom catalysts that sit below the dashed line are selective for the CO2RR.

Reaction Mechanisms

In this section, key intermediates for the two‐electron transfer products in the CO2RR are presented. We assume that, based on BEP relations, lower reaction energies lead to lower energetic barriers and, therefore, the relative reaction free energy of intermediates are able to determine the selectivity of a material toward a product.[ , ] Here, the free energy diagram of Al (Figure 2) and relative energy of reaction intermediates (Table S20) are examined in detail while other SACs considered in this study are discussed generally in the main text. Relative free energy of intermediates and free energy diagrams for all systems can be found in the Supplementary Information (Supplementary Note S5, Tables S21–S26 and Figures S3–S8).

Figure 2

Relative free energy diagram for Al‐doped into the a(i) C3‐site and b(i) N4‐site, with respect to *+CO2. * denotes the clean catalyst surface, while the nomenclature, *CO2, for example, indicates adsorbed CO2 on the surface. Dashed lines show the *CO formation pathway while solid lines show the *HCOOH formation pathway. Inset is a top‐view of the site geometry. Lowest energy adsorbate binding configurations are shown in a(ii)–(vi) and b(ii)–(vi). Al is shown in teal, C in grey, O in red, N in blue and H in pink.

Al

Al has not been widely studied for catalytic performance. Historically, as a foil, Al was found to be inactive for the CO2RR and highly selective for the HER. However, recently it was reported that Al confined in a metal‐organic framework can both suppress the HER and enhance the selectivity for the CO2RR.

The free energy diagrams for Al/C3 and Al/N4 are shown in Figure 2a(i) and b(i), respectively and the geometries of the stable intermediates found in Figures 2a(ii)‐a(vi) and 2b(ii)‐b(vi). All energies shown are free energies comprising both the binding of the indicated molecule itself (relative to standard atomic references) and the binding of the molecule to the catalytic site. These are then referenced to the unbound CO2 molecule to produce the relative free energy shown (see Supplementary Note S5). In the C3‐site, the initial adsorption of the CO2 molecule is isoenergetic with CO2 with a free energy change of 0.01 eV. From *CO2, two possible reaction pathways exist based on the abstraction of a proton from water, either to form *OCHO or *COOH. Figure 2a(i) shows that the hydrogenation of *CO2 to form *OCHO is stabilised by more than 1.14 eV over *COOH. From the *OCHO intermediate, the addition of a second proton to form *HCOOH is essentially isoenergetic (‐0.03 eV free energy change from *OCHO). In contrast, the *COOH intermediate requires a free energy change of 0.34 eV to take up a second proton and produce *CO+H2O. Finally, the free energy of desorption must be overcome to release the two‐electron transfer products. Here, *CO desorption has a reaction free energy of 0.32 eV while *HCOOH desorption has a reaction free energy of 1.53 eV.

In the Al/C3‐site, the formation of *HCOOH has a U of 0.03 V, with the potential limiting step being the second hydrogenation from *OCHO to *HCOOH. In the reaction mechanism toward *CO+H2O production, the potential limiting step is also the second hydrogenation step, from *COOH to *CO+H2O, with a free energy change of 0.34 eV, resulting in a U of −0.34 V. The strong thermodynamic preference for *OCHO over *COOH suggests that *HCOOH will be the thermodynamically favoured product. The desorption of *CO and *HCOOH are not electrochemical steps, therefore these steps do not determine U.

On the Al/N4 SAC, the adsorption of CO2 is also almost isoenergetic with *+CO2, with a free energy change of 0.12 eV, however as this is not an electrochemical step, it cannot be overcome by applying a potential. However, it is generally accepted that a free energy change of up to 0.75 eV is surmountable under reaction conditions. Therefore, it is feasible that this small uphill step can be overcome. After the reaction of *CO2 with H++e−, the *OCHO intermediate is stabilised by 1.24 eV over *COOH. The second hydrogenation step from *OCHO to *HCOOH and the step from *COOH to *CO+H2O both have a reaction free energy of 0.64 eV. Finally, the energy of desorption of *CO is 0.08 eV while for *HCOOH it is considerably larger at 1.02 eV.

In the N4‐site, for both *CO+H2O and *HCOOH production, the second hydrogenation step is responsible for determining the U, from *COOH to *CO+H2O and *OCHO to *HCOOH, respectively. Both reaction pathways also have the same U of −0.64 V. However, despite the equivalent U for the formation of both *HCOOH and *CO+H2O products, the considerable thermodynamic favourability of *OCHO over *COOH, makes it likely that the pathway to *HCOOH will be thermodynamically favoured.

In both the C3‐ and N4‐site, those intermediates that bind through O (*OCHO, Figures 2a(v) and 2b(v) and *HCOOH, Figures 2a(vi) and 2b(vi)) are considerably more stable than their alternative reaction intermediates, *COOH and *CO, respectively, which bind through C (Figures 2a(iii), 2a(iv), 2b(iii), and 2b(iv)). The preference of p‐block metals for binding O over binding C and, as a result, the preferential binding of *OCHO over *COOH, is a well‐known phenomenon.[ , , , ] The oxophilicity of the p‐block elements is typically attributed to more effective hybridisation between the valence orbitals of the p‐block elements and O compared to C. Here, the difference in binding of the *OCHO and *COOH intermediates on the Al/C3 SAC were examined through projected density of states (PDOS) (Figures S9–S11), electron localisation function analysis (Figure S12) and charge density distribution (Figure S13). The PDOS, in particular, displays considerably better overlap of the Al p‐orbitals with the O p‐orbitals of *OCHO, compared to the C p‐orbitals of *COOH. In addition to the PDOS, the increased stabilisation of *OCHO over the *COOH intermediate may be due to the *OCHO interacting with the Al atom through two O atoms (Figure 2a(v)) whereas the *COOH intermediate interacts with the Al atom through a single C atom (Figure 2a(iii)). However, we expect that the bi‐dentate‐like binding of the *OCHO intermediate is less influential in the binding strength than orbital overlap. We reach this conclusion because on the Al/N4 system, the *OCHO intermediate is bound through a single O (Figure 2b(v)) but ΔG(*OCHO) is of a comparable magnitude between the C3‐ and N4‐sites.

Other Metals

The free energy diagrams for all of the remaining metals and defect sites considered in this study (Tables S21–S26 and Figures S3–S8) stabilise the formation of *OCHO over the formation of *COOH. In the C3‐site, the *OCHO intermediate is stabilised by a minimum of 0.45 eV when doped with Ga, and a maximum of 1.14 eV when doped with Al, compared to the *COOH intermediate on the same SAC. In the N4‐site, the stabilisation of *OCHO over *COOH is greater than in the C3‐site, being between 0.74 eV when doped with Ga and 1.36 eV when doped with Bi. Thus, we assume that *HCOOH will be the major thermodynamic product and the predominant mechanistic pathway for these SACs. Therefore, only results toward *HCOOH formation are discussed in the remaining sections, with the pathway toward CO formation being shown in the Supplementary Information (Figures S3–S8). In addition, we highlight that with only one exception (Au/C3), the U is lower for the pathway toward *HCOOH than for the pathway toward *CO+H2O, indicating the reaction toward *HCOOH formation will be initiated at a lower potential compared to the *CO+H2O product. Our finding, that p‐block doped SACs likely proceed via the *OCHO intermediate toward *HCOOH formation is consistent with findings reported for p‐block metal surfaces both experimentally[ , , , , ] and computationally.[ , , , ]

Scaling Relations

Scaling relations are integral to computational electrocatalysis, as they provide a predictor of how an intermediate will bind to a catalyst. Scaling relations correlate the binding energy of, for example, a C‐binding intermediate such as *COOH, with a simple C‐bound species such as *CO and have been widely applied to transition‐metal surfaces. However, for transition metal doped graphene‐based SACs, scaling relations have been met with varying success with Yuan et al. finding scaling relations to hold but with Back et al., Siahrostami et al. and Cui et al. finding that scaling relations are broken for their systems. Cui et al. went further to propose that the broken scaling relations resulted from differences in the geometry of the adsorbed species leading to changes in orbital overlap between the adsorbate and the active site. In support of this proposal by Cui et al., similarity in adsorption geometry has been cited elsewhere as the cause of scaling relations being held.

Here, scaling relations for intermediate binding on the SACs considered in this study are explored. The free energy of O‐binding intermediates, ΔG(*CO2) and ΔG(*HCOOH), were plotted with respect to ΔG(*OCHO) (Figure 3). We do not plot ΔG(*OCHO) with reference to itself. The R2 values for these relations are low, being 0.34 and 0.52 for ΔG(*CO2) and ΔG(*HCOOH), respectively. Despite the scaling relations plotted with reference to ΔG(*OCHO) having lower R2 values than when plotted in reference to ΔG(*OH) (Figure S14), ΔG(*OCHO) provides a sensible reference for constructing a limiting potential volcano, because in the two‐electron transfer process toward *HCOOH formation, U is directly related to the difference in binding between ΔG(*CO2) and ΔG(*HCOOH), with respect to ΔG(*OCHO).

Figure 3

Scaling relations of (a) ΔG(*CO2) and (b) ΔG(*HCOOH) with respect to ΔG(*OCHO).

Scaling relations for C‐binding species in the pathway toward *CO+H2O production were also considered with reference to both ΔG(*CO) (Figure S15) and ΔG(*COOH) (Figure S16), however neither reference provided an R2 of >0.38.

Limiting Potential Volcano

From Figure 3, it is clear that the doped SACs do not adhere well to the scaling relations. Despite the poor R2 value of the scaling relations, these relations were used to construct a limiting potential volcano (Figure 4). The methodology for constructing volcano plots is explained in detail by Casey‐Stevens et al. Table 1 shows that most commonly, U results from the second hydrogenation step from *OCHO to *HCOOH. It is interesting to note that those systems that have a ΔG(*OCHO) >0.63 eV are generally limited by the first hydrogenation step from *CO2 to *OCHO (Table 1). One exception to this generality is Pb/C3, which is limited by the second hydrogenation step from *OCHO to *HCOOH (Table 1), despite having a ΔG(*OCHO) of 0.67 eV. Those SACs that have the lowest overpotential (i. e. the difference between U and the equilibrium potential ) are likely the best catalysts for the CO2RR toward *HCOOH. Thus, we can readily identify Sn/C3, Pb/C3 and Pb/N4 as being the best candidate SACs toward *HCOOH formation. A limiting potential volcano was also constructed for *HCOOH production using ▵G(*OH) as a reference (Figure S17) which also identifies these catalysts as promising materials. Furthermore, experimental synthesis of single‐atom Sn doped materials has recently been achieved thus providing an avenue for experimental testing of this computational result. Pb‐doped graphene‐based SACs are yet to be realised experimentally, to the best of the authors′ knowledge.

Figure 4

Limiting potential volcano plot constructed from the scaling relations calculated in Figure 3 using ΔG(*OCHO) as the descriptor. In general, single atom catalysts with a ΔG(*OCHO) <0.63 eV (denoted as the green background) have the second hydrogenation step from *OCHO to *HCOOH as the limiting step, while those single‐atom catalysts with a ΔG(*OCHO) >0.63 eV (denoted as the pink background) have the first hydrogenation step, from *CO2 to *OCHO, as limiting. The exception to this is Pb/C3 which has ΔG(*OCHO) >0.63 eV but has the second hydrogenation step as limiting. The equilibrium potential for the CO2RR toward HCOOH is −0.17 V vs RHE.

Table 1

The theoretical limiting potential (U), overpotential and associated limiting step for the CO2RR toward *HCOOH. The equilibrium potential for this reaction is −0.17 V vs. RHE.

Metal	Site	PDS	UL [V]	Overpotential [V]
Al	C3	*OCHO→*HCOOH	0.03	−0.20
	N4	*OCHO→*HCOOH	−0.64	0.47
Au	C3	*OCHO→*HCOOH	−0.47	0.30
	N4	*CO2→*OCHO	−0.95	0.45
Bi	C3	*CO2→*OCHO	−0.74	0.57
	N4	*OCHO→*HCOOH	−0.35	0.78
Ga	C3	*OCHO→*HCOOH	0.37	−0.54
	N4	*OCHO→*HCOOH	−0.62	0.33
In	C3	*OCHO→*HCOOH	−0.39	0.22
	N4	*OCHO→*HCOOH	−0.50	0.18
Pb	C3	*OCHO→*HCOOH	−0.21	0.04
	N4	*CO2→*OCHO	−0.22	0.05
Sn	C3	*CO2→*OCHO	−0.24	0.07
	N4	*CO2→*OCHO	−0.60	0.43

U [V]

In previous studies of transition metal doped graphene‐N systems, the activity of CO formation was attributed to the nature of the metal dopant in the site, with the defect site being less relevant to catalytic performance. Here, however, we observe that some metals, i. e. Pb and In, have very similar behaviour regardless of the defect site that the dopant sits in, while other metals, i. e. Al and Ga, have highly site‐dependent behaviour.

Activity Descriptors

Scaling relations are an effective method of describing likely catalytic behaviour but do not provide an explanation of why one catalyst outperforms another. Therefore, activity descriptors can be utilised to describe the intrinsic electronic and/or geometric structure of the catalyst to allow prediction of catalytic ability. To date, the development of activity descriptors has solely focussed on describing trends across the transition metals. Here, we examine the applicability of these activity descriptors for p‐block doped SACs.

Unlike transition metal elements, it is very difficult to establish trends across the p‐block elements. However, in theory, the same principles should apply to this different block of elements. That is, the bonding interaction between the SAC and adsorbate should be predictable from some property of the SAC, whether it be electronic structure, geometric structure or a combination of the two. It is well known that the d‐band centre acts as a good proxy for the strength of binding of key intermediates on transition metals.[ , , , , ] However, clearly, the d‐band centre is not an appropriate descriptor for the p‐block elements, which have valence p‐orbitals. Thus, here, we determine the p‐band centre of the optimised clean SACs and examine correlation to ΔG(*OCHO). While we consider the p‐band centre as a descriptor, Au is excluded as it is a d‐block element. The p‐band centre was calculated as:

where ϵ is the energy of the p‐states and are the number of p‐states at energy ϵ. Here, the average p‐band centre was calculated between −10 eV and 5 eV, with respect to the Fermi level.

In Figure 5, the relationship between the total average p‐band centre and ΔG(*OCHO) is examined. ΔG(*OCHO) was selected as it is the key intermediate in the formation of *HCOOH, for which p‐block elements are selective. If the N4‐sites are considered, an R2 value of 0.38 is obtained, however if In/N4 is excluded, a considerably higher R2 value of 0.98 is obtained. In the C3‐site, there are no clear outliers, and the correlation between the average p‐band centre and ΔG(*OCHO) has an R2 of 0.002. To further increase the resolution of this plot, the p‐, p‐ and p‐band centres were calculated in the same way as the average p‐band centre (equation 12) and plotted with respect to ΔG(*OCHO) (Figures S18–S20), with R2 being no greater than 0.20 for any descriptor. The p‐band centre was also calculated over all p‐states and also only below the Fermi level however, no improvement of the relation was observed. As such, it is clear that the p‐band centre is not a robust descriptor for binding energies of intermediates on p‐block doped SACs. It is not altogether surprising that the p‐band centre is not as effective as the d‐band centre because the p‐bands are known to be more broad than the valence d‐band of the transition metals.

Figure 5

Scaling of ΔG(*OCHO) with the average p‐band centre of the dopant element in single‐atom catalysts.

After the inability of the p‐band centre to describe the binding energies of ΔG(*OCHO), previously developed universal descriptors for the catalytic potential of transition‐metal doped C‐based SACs were applied to the present systems to determine if these descriptors could be effective at describing the behaviour of p‐block element doped SACs. Specifically, we examined the descriptors developed by Gong et al. (Figure S21), Xu et al. (Figure S22) and Guan et al. (Figure S23). None of these descriptors provided an R2 above 0.18.

Clearly, the p‐block doped systems exhibit considerably less predictable behaviour than their transition‐metal counterparts. Not least among the issues of predicting the behaviour of the p‐block elements by applying universal descriptors is the issue of valency. The p‐block elements are effected by the so‐called ‘inert‐pair’ effect,[ , ] whereby the difference in energy of the valence s and p orbitals results in the s electrons not being involved in binding for the lower periods (period >4), while in the upper periods (period <4), both s and p electrons are available for bonding. Ga sits on the cusp of the transition from bonding to non‐bonding valence s 2 electrons and, therefore, may exhibit a transitional behaviour between the two states.

Alternative Products

While only the two electron transfer products are considered in this study, the CO2RR is a notoriously complex reaction pathway with many possible electrochemical products. The *OCHO intermediate can react to form *OCH2O, while the *HCOOH intermediate can also react further to form *CHO. Both *OCH2O and *CHO are key intermediates in the reaction pathways to form H2O, CH3OH, CH2O, and CH4.[ , ] In Figure 6, the free energy required to cleave *HCOOH from the surface is displayed. The black dashed line indicates the 0.75 eV which is generally considered to be the energetic change that is surmountable under reaction conditions for non‐electrochemical steps. Systems that sit below this line are able to desorb *HCOOH as the product because the free energy required for this to occur is attainable under reaction conditions. In contrast, systems that sit above this line indicate that *HCOOH is bound sufficiently strongly and, therefore, may react further to form higher order products.[ , ] Those systems that bind *HCOOH sufficiently strongly to react further are Al/C3, Al/N4, Au/C3 and Ga/C3.

Figure 6

The energy of thermal desorption of *HCOOH from the single‐atom catalyst. The black dashed line shows the 0.75 eV threshold that is generally considered to be a surmountable thermal barrier under reaction conditions.

As a final note, we highlight that in this study only the formation of thermodynamic products are considered. Without the explicit calculation of transition state energies and microkinetic modelling, conclusions about the rate at which products form are unable to be drawn. The kinetic product formation is also a vital contributor to catalyst performance and should be considered in future studies.

Conclusions

In conclusion, we carry out a screening study of p‐block doped SACs for the CO2RR. For all of the SACs considered in this study, *OCHO is considerably more stable than the *COOH intermediate. Pb/C3, Pb/N4 and Sn/C3 are identified as having the lowest overpotential of the SACs considered for the formation of *HCOOH. Interestingly, some SACs display highly site‐dependent behaviour (e. g. Al and Ga) while other SACs behave essentially identically, regardless of site (e. g. Pb and Sn). Furthermore, the thermal barrier required to desorb *HCOOH is likely too great on Al/C3, Al/N4, Au/C3 and Ga/C3 to be overcome under reaction conditions, making these catalysts candidate materials for producing alternative products such as CH4 or CH3OH. However, more complete mechanistic studies need to be carried out on these materials to determine what these products might be. Finally, we discover that the p‐band centre is not an accurate descriptor for ΔG of key intermediates, nor are previously developed universal descriptors. As such, more work is required to develop effective descriptors for p‐block doped SACs to describe trends in their catalytic performance.

Conflict of interest

The authors declare no conflict of interest.

As a service to our authors and readers, this journal provides supporting information supplied by the authors. Such materials are peer reviewed and may be re‐organized for online delivery, but are not copy‐edited or typeset. Technical support issues arising from supporting information (other than missing files) should be addressed to the authors.

Supporting Information

Click here for additional data file.