Computational Screening of Doped Graphene Electrodes for Alkaline CO2 Reduction

The electrocatalytic CO2 reduction reaction (CO2RR) is considered as one of the most promising approaches to synthesizing carbonaceous fuels and chemicals without utilizing fossil resources. However, current technologies are still in the early phase focusing primarily on identifying optimal electrode materials and reaction conditions. Doped graphene-based materials are among the best CO2RR electrocatalysts and in the present work we have performed a computational screening study to identify suitable graphene catalysts for CO2RR to CO under alkaline conditions. Several types of modified-graphene frameworks doped with metallic and non-metallic elements were considered. After establishing thermodynamically stable electrodes, the electrochemical CO2RR to CO is studied in the alkaline media. Both concerted proton-coupled electron transfer (PCET) and decoupled proton and electron transfer (ETPT) mechanisms were considered by developing and using a generalization of the computational hydrogen electrode approach. It is established that the CO2 electrosorption and associated charge transfer along the ETPT pathway are of utmost importance and significantly impact the electrochemical thermodynamics of CO2RR. Our study suggests an exceptional performance of metal-doped nitrogen-coordinated graphene electrodes, especially 3N-coordinated graphene electrodes.


INTRODUCTION
The extensive use of fossil resources has escalated the emission of green house gases, particularly CO 2 , and disrupted the atmospheric carbon balance. An appealing approach for maintaining this balance is to utilize renewable energy resources but their intermittent nature presents a serious obstacle in the energy conversion and storage applications (Vasileff et al., 2017;Jia C. et al., 2019). In this regard, converting renewable electrical energy into chemical energy through the electrochemical CO 2 reduction reaction (CO 2 RR) has been identified as an efficient solution (Tian et al., 2018;Jia C. et al., 2019). Recently, an integrated electrocatalytic CO 2 RR process has drawn appreciable attention due to its extra-ordinary features in enabling atmospheric sequestration of CO 2 followed by ambient CO 2 RR to synthesize chemicals or fuels (MacDowell et al., 2010;Vasileff et al., 2017;Jia C. et al., 2019). Further flexibility is obtained by using the electrode potential and reaction conditions (pH, electrolyte) to control reaction thermodynamics and kinetics as well as activity and selectivity.
However, achieving efficient CO 2 RR is a challenging task and requires optimization of the electrode material as well as the reaction conditions. There are multiple reasons for this. Firstly, CO 2 is a highly stable molecule as reflected by its high negative reduction potential (−1.90 V vs. SHE) required to drive the first-electron transfer process. To circumvent this, active electrocatalysts need to be developed, where CO 2 can be converted to several different products through sequences of complex, multistep proton-coupled electron transfer (PCET) steps. In addition, as many CO 2 derived products have similar thermodynamic stability, the developed electrocatalyst has to be selective. Controlling selectivity while retaining high activity is the primary goal in CO 2 RR electrocatalysis and requires exquisite control over the complex PCET chemistry, which depends sensitively on the electrode material, electrode potential, pH, and electrolyte. Therefore, optimizing the electrocatalytic performance presents a serious challenge to CO 2 RR-based applications (Schneider et al., 2012;Jia C. et al., 2019).
In the search for ideal CO 2 RR catalysts, numerous metallic electrodes have been thoroughly explored experimentally and computationally (Back et al., 2015a;Back et al., 2015b;Mistry et al., 2014;Peterson et al., 2010;Peterson and Nørskov, 2012;Shi et al., 2014;Zhu et al., 2014;Hori, 2008;DeWulf et al., 1989;Kim et al., 1988;Hori et al., 1986;Hori et al., 1989;Lu et al., 2014;Chen et al., 2012;Gattrell et al., 2006;Akhade et al., 2014;Bagger et al., 2019). However, most of them suffer from one or several of the following shortcomings: low faradic efficiency and selectivity, high overpotential, CO poisoning, high cost, and/or low abundance. The extended metal surfaces are also subject to the d-band center theory and intrinsic thermodynamic scaling relationships between the reaction energies of CO 2 RR intermediates. These features together with the Sabatier principle set constraints on the achievable electrocatalytic performance of metallic CO 2 RR catalysts and hampers the search for efficient electrode materials (Peterson and Nørskov, 2012;Hansen et al., 2013;Back et al., 2015b;Li et al., 2015).
To circumvent these inherent limitations of current CO 2 RR electrocatalysts, several unorthodox or innovative CO 2 RR electrocatalysts have been suggested recently (Qu et al., 2010;Wang et al., 2015;Lu and Jiao, 2016;Sun et al., 2017;Kibria et al., 2019;Nitopi et al., 2019). Carbon-based electrocatalysts and in particular doped graphene sheets, are among the most promising materials and have been widely investigated as potential CO 2 RR electrodes. Graphene electrocatalysts have several attractive physical properties such as high surface area, high electron mobility, high thermal conductivity, high Young's modulus etc. Furthermore, their geometrical flexibility and electronic structure have been suggested to enable escaping the scaling relationships Li et al., 2015). While pure graphene can facilitate outer-sphere electron transfer reactions, edges or other defects are usually required for appreciable electrocatalytic activity (Brownson et al., 2012). Even higher activity and selectivity can be achieved by doping pristine or defected graphene with metal or non-metallic elements (Varela et al., 2018;Jia M. et al., 2019;Wu et al., 2019).
Given the promising CO 2 RR performance of such graphenebased materials (Varela et al., 2018;Jia M. et al., 2019;Wu et al., 2019), several experimental and computational investigations have explored the origin of their electrocatalytic behavior and to search for new electrocatalysts. Most computational studies have focused on identifying the optimal reaction energies of the electrochemical PCET steps but the piling evidence (Ringe et al., 2019;Lee et al., 2020;Vijay et al., 2020) indicates that the non-PCET CO 2 adsorption step is limiting CO 2 RR. The reason for calculations focusing on electrochemical PCET steps is that the computational hydrogen electrode (CHE) model (Nørskov et al., 2004), in its original form can only account for PCET steps. The importance of non-PCET steps, such as CO 2 adsorption, can be highlighted by comparing two recent CO 2 RR studies Guo et al., 2020). Considering only the PCET steps, led to the conclusion that CO 2 RR on a Fe_4N is thermodynamically facile and a potential as small as ∼0.1 V vs. RHE is sufficient to produce COOH (Guo et al., 2020). However, accounting for non-PCET steps led to a very different conclusion as the CO 2 adsorption itself is the potential limiting step that has a high thermodynamic barrier of ∼0.9 eV .
In addition to considering the electrode material, it should be recognized that the electrocatalytic CO 2 RR activity and selectivity are inherent properties of electrochemical interfaces. As the interfacial properties depend sensitively on the electrode material, the electrolyte, pH, and electrode potential (Lu and Jiao, 2016;Nitopi et al., 2019), they can be used for controlling the reaction environment and the electrocatalytic CO 2 RR performance (Pérez-Gallent et al., 2017;Aljabour et al., 2018;Gao et al., 2018;Guo et al., 2018;Ringe et al., 2019). Recent studies have demonstrated that the CO 2 adsorption step is sensitive to these properties and it also determines the CO 2 RR activity and selectivity (Ringe et al., 2019;Lee et al., 2020;Vijay et al., 2020). Furthermore, the faradaic efficiency of CO 2 RR is higher in alkaline conditions where the competing hydrogen evolution reaction (HER) is suppressed. Utilizing highly alkaline conditions facilitates obtaining high current densities at lower overpotentials (Gabardo et al., 2018) and the issue of carbonate formation can be circumvented using gas diffusion electrodes (Malkhandi and Yeo, 2019). However, one has to consider both bulk and local pH which may differ as H + /OH − are consumed or formed at the interface resulting in a pH gradient between the electrolyte and electrode interface (Bohra et al., 2019;Varela, 2020). Furthermore, these local pH changes are found to be sensitive to the used electrode potential and current density (Lu et al., 2020).
Given the advantages of pH control and highly alkaline conditions, it is surprising that previous computational CO 2 RR studies have solely focused on the acidic CO 2 RR. A crucial difference between acidic and alkaline CO 2 RR is the proton donor: in acidic conditions, the proton is donated by the hydronium ion or some other acid, whereas, under alkaline conditions, the solvent (water) is more likely the hydrogen donor. This difference can introduce changes to the reaction mechanism: A fully coupled PCET mechanism is preferred in acidic conditions, whereas, under alkaline conditions, a decoupled transfer of an electron and a proton may become the dominant mechanism (Koper, 2013a). In order to understand alkaline CO 2 RR, one has to examine the possibility of decoupled electron transfer/proton transfer (ETPT) steps but this requires going beyond the original CHE (Nørskov et al., 2004) model, which is applicable to coupled PCET reactions only. This limitation can be circumvented using more general approaches such as general grand canonical DFT methods (Mermin, 1965;Melander et al., 2019;Melander, 2020) and the decoupled CHE (Lindgren et al., 2020). While these methods are applicable to ETPT pathways as well, they are currently computationally too demanding for large scale computational screening studies and more tractable methods need to be developed.
In this study, we address the CO 2 RR-to-CO on several graphene electrodes in alkaline conditions and account for both PCET and ETPT pathways. To enable this, the commonly applied computational hydrogen electrode method is extended to study the decoupled ETPT reactions. In particular, we consider the effect of partial ET and the potential-dependency of CO 2 adsorption. On this basis a computational screening approach is developed and applied to identify promising doped graphene electrocatalysts for alkaline CO 2 RR. A four-level selection criteria is introduced to rank and select catalysts according to their thermodynamic stability, ability to bind and activate CO 2 , electrosorption properties, and theoretical limiting and overpotentials. These principles allow us to identify thermodynamically stable doped graphene electrodes with low limiting and overpotentials in alkaline environments. Our study outlines that highly exergonic CO 2 adsorption associated with high degree of electron transfer is detrimental for catalytic performance. We find that N-coordinated Pt and Ru-doped electrodes are promising candidates as pH universal CO 2 RR electrodes.

Modeling of Graphene Sheets
Pristine graphene is modeled using a simplified non-periodic structure with 42 carbon atoms forming a honeycomb structure with 14 benzene rings. The dangling carbon atoms at the edges of the sheet are terminated with hydrogen atoms . The flake-based graphene models have been previously shown to be perform well compared to their respective periodic models (Lazar et al., 2013;Shang et al., 2020). This is attributed to the accounts of non-clustering of carbon atoms and negligible edge effects due to hydrogen-terminations, which make the graphene-flake models trustworthy. In addition, terminating the dangling edge carbon atoms with hydrogen lead to uniform and delocalized distribution of charge and orbitals all over the catalyst surface (Tachikawa and Kawabata, 2011), which prevents any possible artificial localization effects. The grapheneflake models are computationally less expensive compared than extended models while providing a faithful description of graphene. Finally, the graphene-flake model allows us to carry out the large-scale screening studies at the hybrid functional DFT-level (see below for details), which is required to accurately capture the adsorption energies of the CO 2 RR intermediates on doped graphene catalysts (Vijay et al., 2020).
Pristine graphene sheets are chemically inert, but their reactivity can be tuned with dopants and defects. The single vacancy is arguably the simplest point defect in graphene; however, larger point defects such as di-and tri-vacancies are also observed in experiments (Warner et al., 2012;He et al., 2014) forming after the coalescence of two or three neighboring monovacancies (Trevethan et al., 2014). Different vacancy and doping structures can be experimentally realized using, e.g., selective bombardment followed by ion deposition (Wang et al., 2012), or more refined chemical synthesis (Varela et al., 2018). These vacancies or defects can host a variety of dopants (Wang et al., 2012;He et al., 2014), which we modeled by introducing foreign atom(s) into the created vacancy. These atoms include platinum group metals (Rh, Pd, Pt, and Ru), coinage metals (Ag, Au, and Cu), base metals (Al, Fe, Ni, Mo, Co, Mn, Zn, and Cr), and non-metals (B and N). Structural information for all the studied systems is provided as a Supplementary Material set.
We consider three different vacancy structures: single-vacancy (SV), di-vacancy (DV), and tri-vacancy (TV), by removing one, two, and three central carbon atoms, respectively. Three SV systems are investigated namely, M_SV, M 2 _SV, and MPt_SV (see Figures  1A-C). The M_SV structure results from the deposition of a single metal on a single-vacant structure and presents the simplest doped structure. In the case of M 2 _SV and MPt_SV structures, homo-and heteroatom dimers are placed on SV as models for nucleated reaction centers (Ferrante et al., 2016). We also investigated an experimentally inspired M 2 _2SV structure (He et al., 2014), where two neighboring single vacancies are both filled with identical atomic dopants, see Figure 1D. The M_DV structure, shown in Figure 1F has a single dopant in a di-vacancy. The trapping of dopant pairs over tri-vacant surface is energetically more favourable than a single metal dopant and have been observed experimentally (He et al., 2014), and for this reason the M 2 _TV structure is considered as well ( Figure 1H).

Density Functional Theory Methods
All density functional theory (DFT) (Hohenberg and Kohn, 1964;Kohn and Sham, 1965) calculations were carried out using Gaussian 09 (RB.01) (Frisch et al., 2009) package with the help of the Gauss View 5 (Dennington et al., 2009) visualizer. Exchange-correlation effects were described using the B3PW91 (Perdew et al., 1991) functional with Grimme's D3 (Grimme et al., 2010) dispersion corrections. The metallic atoms were treated within the effective core potential (ECP) formalism and the LANL2DZ (Hay and Wadt, 1985a;Hay and Wadt, 1985b;Wadt and Hay, 1985) basis set, while for all the other atoms (C, H, O, N, and B) the 6-31 + g (d,p) (Petersson et al., 1988;Petersson and Al-Laham, 1991) basis set was used. The atomic structures were relaxed until the maximum residual force reached below 0.02 eV Å −1 . The satisfactory performance of the adopted computational approach has been demonstrated previously for graphene  and metallic clusters (Verma and Kishore, 2017;Verma and Kishore, 2019).
Both implicit and explicit solvent effects were mostly excluded from the present study. Considering explicit solvation would have been computationally intractable and the implicit solvent models do not markedly affect the adsorption properties of CO 2 RR intermediates in constant charge calculations as, e.g., explicit hydrogen bonding cannot be captured with such models . The exclusion of solvation effects from modeling the electrocatalytic reactions is a common approximation that has been shown to successfully describe a variety of electrocatalytic systems (Frydendal et al., 2015;Busch et al., 2016;Calle-Vallejo et al., 2017;Valter et al., 2018). However, we have checked the effect of implicit solvation for some of the consider structures. In particular, we addressed implicit solvent effects using the polarizable continuum model (PCM) (Scalmani and Frisch, 2010) for the best performing electrodes identified from the vacuum calculations. For this, we computed the overpotentials using Equation 6 with and without implicit solvent for the best performing catalysts, and the differences are at largest 0.16 V vs. RHE corresponding to 0.16 eV difference in reaction energies as shown in Table 1.
Metal-modified graphene structures are subject to spininversion and, therefore, spin-unrestricted DFT was used and a careful investigation of spin multiplicities was necessary to locate the lowest energy spin states. To ensure the correct spinmultiplicity, single point energies at different spin multiplicities were calculated followed by re-optimizations of the atomic structures at the correct spin-multiplicity. The Mulliken charge (q M ), magnetic moment (m), and the most stable spin multiplicity for each dopant in the considered graphene frameworks can be found in Supplementary Table S1. Note that all the formation and binding free energy calculations were performed using the most stable spin state structures and their energetics.
Vibrational frequency calculations were performed for all the optimized structures to confirm that they are true minima. The vibrational frequencies were computed to determine the vibrational entropy and zero-point energy contributions in the thermodynamics of CO 2 RR. The thermodynamic parameters were calculated at 298.15 K temperature and 1 atm pressure. In the Supplementary Material, we present the methodology utilized to compute the thermodynamic entropy, enthalpy, and free energy. The charge analysis was based on the Mulliken scheme (Mulliken, 1955), which is assumed to be accurate enough for the comparison of similar models and similar basis sets.

Computation of Catalyst Stability
Electrode stability is a prerequisite for maintaining catalytic activity for extended operation time (Krasheninnikov et al., 2009;Wang et al., 2013;Back et al., 2017). In doped graphene electrodes, the deactivation is thought to take place through dopant dissolution (Shao et al., 2019). The dissolution will result in the formation of a graphene vacancy, while the dissolved atoms will likely form metallic nanoparticles or other stable products. Hence, the thermodynamic stability of the studied catalysts was evaluated against pristine graphene and the most stable (bulk) phase of the dopant. The catalyst stability was first referenced against gas-phase metal ions and subsequently against the dopant's most stable bulk phase by utilizing experimental cohesive energies (Kittel, 2004). This way the formation energy is compared against the most stable bulk phase of a given dopant. Finally, the formation free energy (G FE ) is computed as: where N is the number of dopant atoms, G MGr is the free energy of the doped graphene electrode, G M is the free energy of the metal atom in the gas-phase, G vacant Gr is the free energy of the vacant graphene, and G coh is the cohesive energy of the dopant(s) (Kittel, 2004).

Computation of Reaction Free Energies
The CO 2 adsorption free energy (G Ads ) over the doped graphene electrodes was calculated using: where G MGr+CO2 is the free energy of CO 2 adsorbed over modifiedgraphene, G MGr is the free energy of the bare modified-graphene, and G CO 2 is the free energy of gas-phase CO 2 .
To model alkaline conditions, we consider water as the hydrogen donor as commonly assumed for alkaline CO 2 RR (Yin et al., 2019). In general, the reaction may take place either via a sequential ETPT or a concerted PCET mechanism (Dunwell et al., 2018). In the ETPT pathway, the first ET denotes partial charge transfer taking place during CO 2 adsorption on the graphene catalyst, which is followed by further partial ET and PT. For both the sequential and concerted pathways, the elementary steps for CO 2 RR to CO are: where Gr (also referred as "*") is the graphene catalyst and δ denotes the partial charge. When Equations 3a,b, are separated, the mechanism is sequential (ETPT) whereas if they are summed the PCET mechanism is followed.
To compute the electrochemical CO 2 RR thermodynamics along both the ETPT and PCET pathways in alkaline conditions, the CHE method was extended to include partial charge transfer and ETPT steps. Note that we are interested in ET taking place during adsorption, i.e., electrosorption, not outersphere ET to a dissolved CO 2 forming CO − 2 (aq). The crucial difference is that the solution phase ET to CO 2 can be considered as an outer-sphere ET reaction the thermodynamics of which are not affected by the electrode material. In that case, the tabulated (Bratsch, 1989) reduction potentials can simply be used for computing the thermodynamics of this outer-sphere ET reaction. For such outer-sphere reactions, the thermodynamics can also be accurately computed using standard approaches (Marenich et al., 2014). Here, we instead consider an innersphere ET taking place during adsorption, and in this case, the electrode cannot be neglected due to the strong interactions and hybridization between the electronic states of electrode and CO 2 . For such reactions, the potential-dependent thermodynamics require going beyond standard approaches and methods such as grand canonical DFT are required (Hörmann et al., 2020). Grand canonical DFT as implemented presently is, however, too costly for large-scale screening studies and we propose a computationally feasible extension to CHE to enable addressing (partial) ET. The extension is motivated by the recent decoupled hydrogen electrode method (Lindgren et al., 2020), which combines the CHE with grand canonical DFT. We utilize a computationally more feasible approximation based on the electrosorption valency (Schmickler and Guidelli, 2014) to model the partial ET taking place during CO 2 adsorption to capture the potentialdependency of this step. By combining the electrosorption valency and CHE, the reaction thermodynamics for each elementary step can be calculated using: is the water dissociation free energy (0.83 eV at pH 14) (Bratsch, 1989) Equations 4a,b are written for the ETPT pathway and when added together, the reaction free energy for the concerted PCET step is obtained, i.e., ΔG 1+2 (RHE) ΔG 1 (RHE) + ΔG 2 (RHE). c (zΔG ads (U)/zU) is the electrosorption valency (Schmickler and Guidelli, 2014), which changes the adsorption energy as ΔΔG ads (U) cΔU. Note that c ≤ δ and the partial charge transfer upon adsorption do not directly change the Gibbs free energy. As discussed in the SM, γ is a potential-dependent quantity and could be obtained from experiments or using potentialdependent grand canonical DFT, (Melander et al., 2019;Hörmann et al., 2020) but such data is not available for the catalysts studied here. Therefore, we assume that the electrosorption valency is independent of the potential (c(U) ≈ c), which is a good approximation at potentials close to the reference potential but breaks down as U ≫ U ref . This is seen as unphysically large equilibrium potential values for CO 2 adsorption steps-the true upper bound for the CO 2 →CO 2 − should be close to its thermodynamic reduction potential of −1.9 V vs. SHE (Bratsch, 1989). Here, we approximate the electrosorption valency using the fitting approach (Schmickler and Guidelli, 2014), which gives: where χs are Pauling electronegativities. χ Gr ≈ 2.5 for all metaldoped graphene electrodes. χ CO 2 ≈ 1.5 is computed using the experimental (NIST database, 2020) CO 2 electron affinity and ionization energies for computing its Mulliken electronegativity (Mulliken, 1934). The Mulliken electronegativity is then converted to the Pauling scale using the linear correlation (Herrick, 2005) between the two scales.
To compare and analyze the electrochemical thermodynamics, we also define the thermodynamic equilibrium potential (U eq ), the thermodynamic limiting potential (U L ), and overpotential (η) (Durand et al., 2011) as further discussed in the SM. However, our definition will take into account the ETPT steps with partial charge transfer in terms of the electrosorption valency.
where n e,tot 2 is the total number of electrons transferred in the CO 2 RR to CO.

Structures and Stabilities
Thermodynamic stability against dissolution is a key material property mandatory for maintaining electrocatalytic activity and, therefore, we first addressed the stability of our graphene model electrodes against pristine graphene and the thermodynamically most stable (bulk) phase of the dopant. Figure 2 summarizes the Gibbs free energies for formation, computed according to Equation 1. Additional data on formation thermodynamics with numerical values, zeropoint energy corrections, and formation enthalpies are reported in the Supplementary Tables S2-S5. We found that majority of the studied dimers (M 2 _SV and MPt_SV) in single vacancies are thermodynamically unstable and these results are therefore presented and discussed solely in the SM. For other systems, the findings of Figure 2 are discussed in more detail below.

Single and Two Mono-Vacant Graphene
The M_SV graphenes are the simplest doped graphene structures (see Figure 1A). Calculations showed that all metal atoms exhibit slight protrusion out of the graphene plane upon bonding with the dangling carbon atoms, and the vertical distance between a metal atom and graphene plane varies from 0.73 Å to 1.19 Å in good agreement with previous computational studies (Li et al., 2015;Zhao et al., 2019). About half of the M_SV systems are thermodynamically stable as can be seen from Figure 2. The free energy of formation ranges from −0.19 to −6.47 eV; however, doping with coinage metals and a few of base and platinum group metals is thermodynamically unfavorable (see Supplementary Table S5). The formation free energies of Band N-doped graphene are above −6 eV, which demonstrates their extremely strong interaction with the graphene single vacancy.
In the case of M 2 _2SV, the metal dopants are also elevated out of the graphene plane similar to M_SV. Almost all M 2 _2SV systems are symmetric with the dopants occupying two vacancies. The formation free energy of stable M 2 _2SV electrodes ranges from −0.42 eV to −5.70 eV. Doping by nonmetallic B and N atoms is extremely exergonic (∼−5.5 eV). We found that, upon structural optimization, the Ag, Al, and Cu dopants spontaneously drift away from the vacancy to form dimers, therefore, they are excluded from electrochemical thermodynamics studies.

Di-and Tri-vacant Graphene
In the M_DV graphene, dopant atoms are bound to four nearest undercoordinated C atoms and remain in the basal plane due to the large size of the double vacancy compared to the dopant's atomic radius. The considered non-metallic elements B and N feature benzene-like hexagons with two out of the four undercoordinated carbon atoms. According to the formation free energy analysis, all M_DV graphene structures are highly stable with formation energies ranging from −1.43 eV to −5.28 eV.
The tri-vacancy of M 2 _TV models can accommodate the dopant pair in two different orientations: 1) both dopant atoms above the tri-vacancy (Supplementary Figure S4A) and 2) one atom above and one below the tri-vacancy (Supplementary Figure S4B). We found that the latter geometry is more stable (see Supplementary Figure S4C) and was selected for further studies. In the M 2 _TV structure, the dopant atoms are bound to the three nearest undercoordinated carbon atoms. Non-metallic B and N dimers break down during structural optimization to form benzene-like structures, and the resulting structures mimic mono-vacant graphenes. The computed formation free energies show that non-metallic dopants in the M 2 _TV structure are highly stable, while metallic M 2 _TV electrodes are at best modestly stable, and Ag-, Au-, Cr-, Mo-, and Ru-doped graphenes are unstable.

N-Coordinated Graphene
All the optimized M_3N and M_4N structures are highly symmetric and show no buckling behavior. We find that the average M-N bond distances along with other geometrical features of both N-coordinated structures are similar to their carbon coordinated counterparts (M_SV and M_DV) Kattel et al., 2012). Our calculations suggest that the incorporation of three nitrogen atoms enhances the dopant's stability as can be inferred from the formation free energies ranging from −1.22 eV to −3.28 eV, which are more exergonic than their M_SV counterparts. The formation free energies of M_4N structures vary from −1.84 eV to −2.99 eV and are comparable to the values computed for M_DV models. However, it is evident that the non-metallic dopants prefer 4C-coordination as their formation free energies (see Figure 2) are much more exergonic than in the 4Ncoordination environment. Therefore, the 4N-coordination environment does not enhance the stability of metallic dopants compared to 4C-coordination in agreement with previous computational results . In general, the geometric and energetic stabilities of modifiedgraphene sheets depend strongly on the dopant and vacancy but some overall trends are observed. The Ag-, Au-, Cr-, and Modoped carbon-coordinated graphene structures are always unstable irrespective of the type of vacancy whereas all the considered dopants in 3N-and 4N-environments are highly stable, which is in line with previous studies (Tripkovic et

CO 2 Adsorption
As discussed in the Introduction, CO 2 adsorption and activation are crucial steps for CO 2 RR in alkaline conditions as the adsorption is often thermodynamically unfavorable and limits CO 2 RR activity (Jones et al., 2014;Gauthier et al., 2019;Yu et al., 2019;Vijay et al., 2020). Furthermore (partial) charge transfer during CO 2 adsorption may take place and this process is sensitive to the electrode potential but formally independent of solution pH. CO 2 may also be activated upon adsorption, which is seen as the deviation of the O-C-O angle from 180 o . Under alkaline conditions, where only weak proton donors (water) are present, the CO 2 adsorption and an accompanied ET step are important contrary to pure PCET steps. Therefore, we closely monitor the CO 2 adsorption geometry and energy over all of the considered electrodes, regardless of electrode stability. Figure 3 displays CO 2 binding geometries on the stable M_SV, M 2 _2SV, M_3N, M_DV, and M_4N electrodes whereas binding geometries on stable M 2 _SV, MPt_SV, and M 2 _TV structures are illustrated in Supplementary Figure S5. Electronic structure information such as magnetic moments and stable spin multiplicities are supplied in Supplementary Table S6 for each system, and numerical values for adsorption energies, enthalpies, and free energies are provided in the Supplementary Tables S7-S10.
Among the M_SV electrodes, CO 2 adsorption is exergonic only on Fe_SV (−0.12 eV). In this case, CO 2 remains linear upon adsorption (see Figure 3A) and the oxygen atom interacts weakly with the metal site having Fe-O distance 2.16 Å and small electron transfer from CO 2 (0.22 e). A similar pattern is observed for the M_DV electrodes as well. The CO 2 adsorption is exergonic only on Fe_DV (see Figure 3D) with similar binding configuration and free energy as observed for Fe_SV; however, the charge transfer (0.47 e) from CO 2 is rather high compared to Fe_SV. On the B-and N-doped SV sheets, CO 2 merely physisorbs as seen from the long B-C and N-C distances, which are over 3.1 Å. This is similar to the CO 2 adsorption on the pristine graphene (Rad and Foukolaei, 2015). The activation of CO 2 , i.e., non-linear CO 2 structure, occurs only on Mn_SV and Cr-and Mo-based SV and DV structures but the CO 2 adsorption is thermodynamically unfavorable. Note that the formation of Crand Mo-doped SV electrodes are thermodynamically unstable as well. Figure 3B illustrates the adsorption geometry of CO 2 on Mn 2 _2SV. The adsorption geometry is similar on all the CO 2activating M 2 _2SV structures. CO 2 adsorption is exergonic on B-, Co-, Mn-, and N-based 2SV structures with adsorption free energies varying from −0.18 eV to −0.68 eV. However, only Co 2 _2SV and Mn 2 _2SV are able to activate CO 2 and facilitate modest electron transfer to the molecule. While Fe 2 _2SV and Ru 2 _2SV structures can also activate CO 2 , they exhibit thermoneutral adsorption. On the other hand, upon structural optimization, the coinage and a few base metal (Al, Ni, and Zn) dopants stem out of the vacancies and form CO 2 -dimer complex above graphene structure, and hence these systems were excluded from further studies.
Apart from Ag-and Cu-doped M_3N models, the majority of M_3N structures exhibit promising CO 2 adsorption characteristics as CO 2 both adsorbs and becomes activated with the O-C-O angle varying between 127.9 o to 150.4 o . The molecule interacts with the metal dopant via the carbon atom and one oxygen in a bidentate conformation shown in Figure 3C for the Mo_3N structure. In the case of non-metallic dopants, CO 2 favors a monodentate adsorption configuration. The free energies of adsorption range from −0.24 eV to −0.74 eV being comparable to the values found for the M 2 _2SV structures. CO 2 adsorption on 4N-coordinated electrodes is surprisingly unfavorable: roughly half of the M_4N structures can activate CO 2 but only Mo_4N (see Figure 3E) exhibits exergonic adsorption free energy (−0.89 eV) with exceptionally large charge transfer (−0.67 e). Similar to M_3N graphenes, CO 2 favors a bidentate adsorption geometry on Mo_4N with O-C-O angle of 130.7 o . Given the experimentally proven performance of 4N catalysts, it is unexpected to observe rather high endergonic adsorption energies on these electrodes. Possible reasons for the observed weak adsorption include, e.g., the neglected axial ligands and dipolar interactions as further elucidated in the Discussion section, or the deformation of the electrode structure upon adsorption. To test the latter, we computed the deformation energy for Fe_4N in the same way as in ref (Nykänen and Honkala, 2011). but the contribution from structural deformation is very small and therefore unfavourable

Selection Criteria
After analyzing the electrode stabilities and CO 2 adsorption, we develop the selection criteria for identifying the most promising modified graphene CO 2 RR electrocatalysts from a large pool of materials. For instance, the doping of heteroatom to the graphene electrode should be highly stable to avoid any possible dissolution. The binding of CO 2 to the electrode material should be strong enough to avoid being the limiting step at industrially relevant current densities/potentials. Furthermore, at high pH conditions, Tafel slopes approach 120 mV/dec (Gabardo et al., 2018), which is indicative of electron-transfer (ET) step being the limiting step (Dunwell et al., 2018). Therefore, CO 2 should be activated upon adsorption, become negatively charged, and bind strongly: these features are ranked based on O-C-O angle, negative charge transfer, and negative G Ads values.
To rank and identify prospective CO 2 RR electrocatalysts, we introduce a four-level selection criteria (see the flow chart in Supplementary Figure S6). First, we screen all bare electrodes and select the ones with negative formation free energy but disregard the systems with geometric instability. The electrodes fulfilling the first criterion are subjected to CO 2 adsorption energy screening. Catalysts, over which the adsorption free energy is below or equal to 0.15 eV (i.e., more adsorbing) are selected. Third, after passing the screening of binding free energies, the next step is CO 2 activation for which we set the threshold value of 150 o for O-C-O angle. Fourth, after previous three criteria, we consider the net charge transfer due to the binding of CO 2 to the electrode as another criteria. In this section, the systems demonstrating net charge transfer below or equal to 0.15 are selected. The upper range of 0.15 eV in the binding free energies of CO 2 and in the net charge transfer is considered to account for any possible number sensitivity due to the DFT methods . This is an acceptable condition considering the fact that the electronic structure results often deviate by 0.1 eV due to different DFT functionals or basis sets.
A combined plot of three key selection parameters, i.e., adsorption free energy (G Ads ), O-C-O angle, and charge transfer (q) for each of the electrode categories, namely, M_SV, M_DV, M 2 _2SV, M 2 _TV, M_3N, and M_4N, is supplied in Supplementary Figure S7. Based on the selection criteria, the qualified systems to study the electrochemical thermodynamics are: M 2 _2SV (Co, Mn, Rh, Fe, and Ru), M 2 _TV (Al, Cu, and Zn), M_3N (Al,N,Rh,Cr,Fe,Mo,Pd,Pt,and Ru), and M_4N (Mo and Ru). Thus, a total of 19 electrodes are selected for further study. Note that none of M_SV and M_DV electrodes qualify the present selection criteria.

COOH and CO Binding: Scaling Relations
CO and COOH adsorption energies typically exhibit scaling relations and are therefore often taken as CO 2 RR activity and selectivity descriptors (Peterson and Nørskov, 2012;Koper, 2013a;Kuhl et al., 2014;Nitopi et al., 2019). However, there are indications Li et al., 2015;Back et al., 2017) that the CO/COOH scaling relations and descriptors are not valid for graphene-based electrodes. To study whether scaling relations can be established among the materials studied here, we attempt to find correlations between the CO 2 , COOH, and CO adsorption free energies. This analysis is carried out for the M_3N and M 2 _2SV electrodes only as M_4N and M 2 _TV electrodes contain too few data points. The binding energies are presented in the Supplementary Table S11 and are utilized in the next section to study the electrochemical thermodynamics.
From Supplementary Figures S8, S9, it is evident that a poor correlation exists between the binding free energies of COOH and CO 2 over M_3N (R 2 0.07) and M 2 _2SV (R 2 0.24) electrodes. Similar results are observed for COOH/CO scaling relationships on M_3N (R 2 0.29) and M 2 _2SV (R 2 0.08) electrodes as well. These poor scaling relations are in agreement with previous studies on graphene-based electrodes (Li et al., 2015;Back et al., 2017;Li et al., 2019;Guo et al., 2020), which indicates that graphene electrodes may break scaling relationships making these materials an interesting class of CO 2 RR electrocatalysts. However, it needs to be recognized that the number of materials subjected to the scaling relation studies is small and too far-reaching conclusions on breaking scaling should be avoided. Yet, it has been shown that even slight modifications to N-doped graphene electrodes lead to different scaling relations (Guo et al., 2020).

Electrochemical Thermodynamics
To address the electrochemical performance of the qualified electrodes, electrochemical thermodynamics for CO 2 reduction to CO were computed using Equations 4, 6. While Figure 4 displays the thermodynamic profiles for selected materials with varying behavior along the PCET and ETPT pathways, the limiting potentials and overpotentials are illustrated in Figure 5. For other materials, not shown here, the potential energy surfaces are displayed in Supplementary Figure S10 and the explicit potential values are collected in Supplementary Table S12.
As displayed in Figure 4A, a large potential difference (0.69 V) between the ETPT and PCET pathways is observed for the Pt_3N catalyst. We attribute this difference to the inclusion of the CO 2 binding in the former mechanism. The CO 2 adsorption step itself is rather exergonic and associated with a modest charge transfer (∼−0.3 e) from the electrode. Along the PCET pathway, COOH and CO formation steps are mildly exergonic at zero potential resulting in positive equilibrium potential. This observation suggests that PCET pathway on Pt_3N is nearly thermodynamically ideal and requires zero limiting potential (blue PES). On the other hand, the ETPT mechanism is far from ideal and a large negative limiting potential is required due to the highly endergonic CO 2 →COOH step and the significant partial charge transfer (∼0.7 e) in this step.
CO 2 RR thermodynamics on the Al_3N catalysts, shown in Figure 4B, provides an example where both the ETPT and PCET pathways attain same limiting potential. Herein, CO 2 adsorption energy is modestly exergonic and features insignificant charge transfer. Furthermore, CO formation is highly unfavourable compared to COOH and is therefore identified as the potentialdetermining step. As the reaction free energy for CO formation (see Equation 4c) is same for both PCET and ETPT pathways, the resulting thermodynamic potentials are also equal.
In Figure 4C, the Ru 2 _2SV electrode is considered. CO 2 adsorption is almost thermoneutral at zero potential (red line) associated with a small charge transfer (∼−0.2 e) from the electrode. The formation of COOH is identified as the limiting step owing to its high endergonicity. Altogether these factors make ETPT and PCET limiting potentials rather similar, which in turn implies that this material is a promising candidate for both ETPT and PCET pathways, i.e., for both alkaline and acid conditions. However, the CO adsorption free energy (∼−1 eV) on the electrode is too strong, which makes it susceptible to catalyzing further reduction or poisoning. The last example, given in Figure 4D, features very different limiting potentials for the ETPT and PCET pathways on the Mo_4N electrode. In this case, CO 2 adsorbs strongly, significant charge transfer (∼−0.7 e) to the molecule takes place, and the COOH formation is quite endergonic at zero potential (red line). Together these factors result in a large negative limiting potential along the ETPT pathway, while the PCET exhibits a fairly modest overpotential. Similar behavior is seen for the Mo_3N electrode as well and we conclude that the N-Mo-modified graphene catalysts are unsuitable for a ETPT mechanism but could work for the PCET pathway.
Comparison of the ETPT potential energy profiles given in Figure 4 shows that CO 2 adsorption and associated charge transfer are pivotal for this mechanism to be operational, but they are absent from PCET due to the assumption of simultaneous electron and proton transfer. Overall, these differences modify the reaction thermodynamics, and we observed that the charge transfer during adsorption has a profound effect on the electrocatalytic thermodynamics as  We found that the slightly exergnonic CO 2 adsorption is beneficial as it enables the formation of a stable CO 2 -catalyst complex and prevents excessively endergonic COOH formation. Furthermore, charge transfer associated with CO 2 adsorption should be as low as possible: if the COOH formation requires only a very small degree of charge transfer, a very high reductive potential is needed to make this step thermodynamically favorable as shown in, e.g., Figure 4D. In a case where no charge transfer takes place during CO − 2 + H + → COOH, the COOH formation step is fully chemical and cannot be controlled by the electrode potential. The very large limiting potentials along the ETPT pathways are likely an artifact of the simple charge transfer and electrosorption model adopted herein, but we believe that the present model enables qualitatively correct comparison of thermodynamics between ETPT and PCET processes.
The thermodynamic potentials are extensively used for ranking or comparing the expected performance of different catalysts, and limiting potentials and overpotentials should be close to zero for the best performing catalysts. Figure 5 presents both potentials for the PCET and ETPT pathways. For all the qualified electrodes, the overpotentials are within 0.7 V (vs. RHE) for PCET pathways. The best-performing electrodes include Ru, Rh, and Mo-doped 3N graphene and Ru and Mo-doped 4N graphene for which the limiting potentials vary from 0 to −0.52 V vs. RHE and overpotentials range from 0.06 to 0.21 V vs. RHE. For the ETPT mechanism, the best-performing electrodes contain Pd_3N, Ru_3N, Pt_3N, Ru 2 _2SV, and Ru_4N with limiting potentials between −0.48 and −0.86 V vs. RHE and overpotentials varying from 0.28 to 0.8 V vs. RHE. The limiting potentials found for the studied structures are either comparable or even smaller than those previously computed for 4N and 3N catalysts (Guo et al., 2020;Pan et al., 2020), non-metallic defected graphene (Siahrostami et al., 2017), and single-atom doped metal catalysts (Lim et al., 2014). Overall, the ETPT mechanism leads to higher overpotentials and limiting potentials than the PCET pathway, which we take to indicate that pure ET steps are detrimental for the reaction thermodynamics as commonly suggested (Koper, 2013b). Interestingly, we find that three of the studied materials-Ru_3N, Pt_3N, and Ru_4N-exhibit very promising electrochemical thermodynamics along both the ETPT and PCET pathways making these materials promising candidates for CO 2 RR under both acidic (PCET) and alkaline (ETPT) conditions. Next, we simulated the three promising ETPT and PCET electrocatalysts (Ru_3N, Pt_3N, and Ru_4N) using the PCM model (Scalmani and Frisch, 2010) for water to check the solvation effects on the overpotentials, and their deviations from the vacuum phase thermodynamics. The adsorption free energies of CO 2 in the aqueous medium are lower (G Ads(aq) Pt_3N −0.49 eV; G Ads(aq) Ru_3N −0.08 eV; and G Ads(aq) Ru_4N 0.02 eV) than those of computed for the vacuum phase, but charge transfer from electrodes to the CO 2 in the solution phase is significantly higher (−0.5 e to −0.6 e). The aqueous phase PCET and ETPT limiting potentials of Pt_3N and Ru_3N are slightly lower compared to vacuum phase, whereas, slightly higher limiting potentials are observed for Ru_4N. As shown in Table 1, including solvation introduces only small changes to overpotentials. The largest change is observed for Ru_4N (Δη 0.16 V), which is still relatively small given typical DFT inaccuracies. We ascribe the observed difference to slightly more stable adsorption of COOH* in the solution phase compared to vacuum phase. Finally, it can be concluded that the implicit solvation framework does not significantly change the thermodynamics computed in gas-phase, and thus the conclusions remain largely unaffected by the implicit solvation.

DISCUSSION
Our computational predictions show that the CO 2 RR activity and selectivity on a given graphene-based electrode depend on the mechanism, i.e., whether the PCET or ETPT pathway is followed. For conditions favoring the PCET mechanism, the N-coordinated graphene electrodes demonstrate the best performance among the studied models. The aforementioned holds true also when ETPT is operational but this typically leads to larger limiting potentials. Pd_3N is an exception to this rule and is the only electrode for which a lower overpotential along the ETPT than the PCET pathway is observed. This anomaly is caused by the slightly positive CO 2 binding energy and minor charge transfer. In general, we find that the coupled PCET is thermodynamically more favorable but for kinetic reasons the ETPT pathway may be preferred (Koper, 2013b).
The obtained results show that the ETPT pathway is highly sensitive to the CO 2 electrosorption energy and charge transfer. Ignoring these important features from mechanistic consideration and focusing solely on the PCET mechanism, often leads to a different potential-determining step, and underestimation of thermodynamic potentials. For example, the current and previous computational studies Vijay et al., 2020) on Fe_4N show that the CO 2 binding itself introduces a thermodynamic barrier of ∼0.9 eV, whereas, the PCET pathway (Guo et al., 2020) presents a quite small (∼0.1 eV) thermodynamic barrier for producing COOH. These results demonstrate that CO 2 adsorption directly affects the CO 2 RR elementary thermodynamics and should therefore be considered. The adsorption process and its potentialdependency are particularly important under alkaline conditions where the coupled PCET becomes less likely, and at high current densities where mass transfer and adsorption are the limiting processes. We identified that slightly exergonic CO 2 adsorption associated with minor charge transfer is a desirable feature for promising catalysts. The conclusions from a large number of studies (see Introduction) neglecting the adsorption step and focusing only on PCET steps are limited to acidic conditions or to the electrodes that cannot catalyze the ETPT pathway.
The importance of the adsorption step can be further illustrated by the widely studied Fe_4N electrode. Our results suggest that the Fe_4N electrode Li et al., 2019;Vijay et al., 2020) does not belong to the best performing CO 2 RR materials. We attribute this to the thermodynamically unfavorable CO 2 adsorption step, which agrees with the recent finding that the adsorption step on Fe_4N is thermodynamically uphill even under high electrode potentials and field strengths (Vijay et al., 2020). Therefore, CO 2 RR on Fe_4N is limited by CO 2 Frontiers in Energy Research | www.frontiersin.org February 2021 | Volume 8 | Article 606742 adsorption and mass transfer at industrially relevant current densities (Vijay et al., 2020). Even though our decoupled pathway and CO 2 RR analyses are based on a relatively simple and ideal electrosorption valency concept, the identification of CO 2 adsorption as a crucial step is expected to be rather general.
In particular, our approach provides a fast approach to identify a catalyst, where ET during adsorption is an important feature. In addition to accounting for ETPT and electrosorption, future studies should also consider the role of "innocent" ligands introduced by the pH or the supporting electrolyte. A recent joint experimental and computational work  has shown that axial ligands are important in determining the CO 2 RR to CO performance on M_4N-type materials. While the experimental results exhibit a high current density and Faradaic efficiency toward CO at applied potentials of −0.5 and −0.6 V vs. RHE, the calculations predict a significant thermodynamic barrier of ∼0.9 eV for CO 2 adsorption in agreement with our results. As the initial computational results were not in line with the experimental findings, the effect of axial adsorption of H 2 O and OH was considered. While the adsorption energy was found to be further destabilized with axial OH adsorption, the axially adsorbed H 2 O stabilized CO 2 adsorption on Fe_4N by ∼0.3 eV. Including the axial H 2 O ligand in the computational Fe_4N model resulted in a better agreement with experiments as the CO 2 adsorption posed only a small thermodynamic barrier (0.06 eV) at −0.6 V vs. RHE. Furthermore, the simulated Pourbaix diagrams  confirmed that Fe_4N electrode may bind axial ligands on both sides. While the kinetic role of the axial ligands was not considered for CO 2 RR, the oxygen reduction reaction kinetics were shown to be very sensitive to the presence of "innocent" axial ligands (Rebarchik et al., 2020).
Previous studies (Peterson and Nørskov, 2012;Kuhl et al., 2014;Nitopi et al., 2019) suggested that CO binding strength is a descriptor that determines the CO 2 RR product distribution. Materials with strong CO adsorption are either poisoned or produce C 1 -or C 2 -species or form hydrogen via the competing HER. On the contrary, electrodes binding CO weakly yield CO as the major product due to favourable CO desorption kinetics. Apart from Al-, Cr-, and Pd-dopants, the majority of 3N-coordinated electrodes bind CO quite strongly suggesting that they might be poisoned by CO. The CO adsorption on Al_3N (η ETPT/PCET 0.69 V) is weak and leads to facile CO desorption. Intermediate CO adsorption energies on Cr_3N and Pd_3N suggest that these two materials are promising electrocatalysts for producing C 1 -or C 2molecules; however, the former suffers from rather high limiting potential making Pd_3N the only 3N-graphene electrode suitable for further CO reduction. With the exception of Ru 2 _2SV, all the other 2SV electrodes exhibit CO binding strength within ∼−0.7 eV indicating that they will primarily produce CO. The M 2 _2SV electrodes show low overpotentials but feature high negative limiting potentials (above −1V) except for Ru 2 _2SV, which is identified as the most promising M 2 _2SV structure. Both M_4N and M 2 _TV electrodes suffer from either strong CO binding or high limiting potentials, and are therefore unsuitable for CO production.

CONCLUSION
In this work, we identified promising modified graphene electrodes for CO 2 RR to CO under alkaline conditions. Robust selection criteria based on thermodynamic stability, CO 2 adsorption thermodynamics, and potential-dependent reaction free energies were devised and applied. The computational hydrogen electrode concept was extended to treat decoupled PCET steps at electrode surfaces to account for the possible decoupled ETPT mechanism in alkaline conditions. We utilized this development to evaluate the effect of pure charge transfer during adsorption, i.e., electrosorption, and found that a high degree of partial charge transfer during CO 2 is detrimental to electrocatalyst performance and that moderately strong CO 2 adsorption energy without charge transfer leads to promising electrode materials. We identified metal sites coordinated to three nitrogen atoms (M_3N) and two single vacancy metal (M 2 _2SV) electrodes as highly promising materials for CO 2 RR following either the coupled or decoupled pathways. N-coordinated Ru and Pt electrodes exhibit promising characteristics for both coupled and decoupled pathways making these materials interesting candidates as pH-universal CO 2 RR electrodes and call for further experimental and computational studies.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be accessed at https:// gitlab.com/ 1341 compcatjyu/co2rr-on-graphene-electrodes.

AUTHOR CONTRIBUTIONS
AV performed the calculations and analysed the data. MM developed the thermodynamic model and oversaw the analysis. KH oversaw the calculations and analysis. All authors contributed to conceptualising the work and writing the manuscript.

ACKNOWLEDGMENTS
The computational resources were provided by the CSC-IT Center for Science, Espoo, Finland (https://www.csc.fi/en/) and the FGCI-Finnish Grid and Cloud Infrastructure.