The geometry optimization and energy calculations were performed within the density functional theory (DFT) framework (Kohn and Sham, 1965) by using the Vienna ab initio simulation package (VASP) (Kresse and Furthmüller, 1996a). The ion–electron interaction was described by the projector-augmented wave (PAW) potentials (Blöchl, 1994). The generalized gradient approximation parameterized by Perdew, Burke, and Ernzerhof was utilized as the exchange-correlation functional (Perdew et al., 1996). The kinetic energy cutoff of 400 eV was adopted for the plane-wave expansion. The armchair-edged ribbon, zigzag-edged ribbon, and periodic graphene slab were sampled with 4 × 1 × 1, 1 × 4×1, and 4 × 4 × 1 Monkhorst−Pack k-point grids (Delley, 2000), respectively. During the geometry optimization, all atoms were relaxed until the total energy was converged to 1.0 × 10^–5 eV/atom, and the force was converged to 0.01 eV/Å. In addition, we considered the van der Waals (vdW) interactions by employing the semiempirical DFT-D2 forcefield approach (Grimme, 2006).

The lattice parameters of 8.52 Å × 24.6 Å and 25.6 Å × 9.84 Å were set to model the armchair-edged graphene nanoribbon (including pyrrolic edge) and zigzag-edged graphene nanoribbon, respectively. The lattice parameters of 9.84 × 9.84Å were adopted to model the periodic graphene slab. Perpendicular to all graphene structures, a vacuum layer of 15 Å was set, which was sufficiently large to minimize the image interactions.

The adsorption energy ( Δ E ads ) of adsorbates was defined as follows: Δ E ads = E substrate+adsorbate − ( E substrate + E adsorbate ) , (1) where E _{substrate+adsorbate} is the total energy of the substrate with adsorbed molecules. E _substrate and E _adsorbate are the energy of the isolated substrate and free molecule, respectively.

The computational hydrogen electrode (CHE) model (Norskov et al., 2004) was adopted to evaluate the free energy change during the CO₂RR process. In the CHE model, the hydrogen atom is in equilibrium with the proton/electron pair at 298.15 K and 1 atm of pressure. In other words, the half chemical potential of gas-phase H₂ is equal to that of a proton/electron pair at 0 V in an aqueous solution.

The Gibbs free energy change (ΔG) for each elementary CO₂RR step involving proton/electron pair transfer was calculated by the expression (Norskov et al., 2004; Zuluaga and Stolbov, 2011): Δ G = Δ E + Δ Z P E − T Δ S + Δ G u + Δ G pH , (2) where ΔE is the change of reaction energy based on DFT calculations. ΔZPE and ΔS are the change of zero-point energy and entropy, respectively. T refers to the temperature (298.15 K). The zero-point energy (ZPE) of adsorbates has been calculated from the vibrational frequencies. For the free molecules (CO₂, CO, HCOOH, CH₄, CH₃OH, etc.) the vibrational frequencies and entropies are obtained from the NIST database (http://webbook.nist.gov/chemistry/). ΔG _U = -neU, where n is the number of transferred electrons, e is the elementary charge of an electron, and U is the electrode potential vs. RHE. ΔG _pH= 2.303 k _BT *pH, k _B is the Boltzmann constant. In this work, the value of pH was set as 0 for the acid medium (Faccio et al., 2010; Shang et al., 2010). Approximate solvation corrections with a dielectric constant of ε = 80 are applied for the simulation of an aqueous environment (Mathew et al., 2019).

In previous reports, the N-doped graphene materials have been widely studied as ORR electrocatalysts, which showed better stability and tolerance to methanol crossover effect than commercial Pt/C catalyst (Geng et al., 2011; Lin et al., 2013; Gong et al., 2009; Qu et al., 2010). Under different temperatures, the synthesizability of each type of the N-doped graphene materials is different. It is relatively easy to synthesize different types of N-doped graphene by controlling the temperature (Lin et al., 2013). The studied structures include five N-doped armchair graphene types, four N-doped zigzag graphene types, in-plane graphitic N (GN), and pyrrolic edge N (PyrroN)-doped graphene. For N-doped armchair graphene, it includes graphitic N (A-GN), pyridinic N (A-pyN), hydrogenated pyridinic N (A-pyN-H), oxidized pyridinic N (A-pyN-O), and pyridinic N hydroxide (A-pyN-OH), as shown in Figure 1. For N-doped zigzag graphene, four structures are considered, i.e., graphitic N (Z-GN), pyridinic N (Z-pyN), hydrogenated pyridinic N (Z-pyN-H), and oxidized pyridinic N (Z-pyN-O). These doped structures could be generated at high temperatures in the pyrolysis process of N-containing compounds (Wu et al., 2011; Li et al., 2012; Wang Q. et al., 2014; Holby et al., 2014).

During the CO₂RR process on the studied compounds, the intermediates mainly include CO₂, COOH, HCOO, HCOOH, CO + H₂O, COHOH, H₂COO, and COH + H₂O. By exploring different adsorption sites (N and its adjacent carbon atoms), the most favorable adsorption configurations and sites are obtained (Supplementary Figures 1–8). Since the two main reactions on various N-doped graphene are HCOOH and CO generation pathways, we focus on the adsorption energies of CO₂, COOH, HCOO, HCOOH, and CO as listed in Table 1, together with the bond distance between the adsorbed intermediates and catalyst surface. It is seen that the adsorption of CO₂ molecule is weak all the time (-0.06 eV∼ -0.13 eV), and linear structure is maintained above the surface. To achieve high selectivity for HCOOH or CO, COOH (or HCOO) should be adsorbed strongly, but HCOOH or CO should be adsorbed weakly for desorption. Therefore, strong COOH (HCOO) binding but weak HCOOH (CO) adsorption is essential for the formation of HCOOH or (CO) (Sharma et al., 2015; Wu et al., 2015).

As shown in Supplementary Figures 2, 3, COOH could not be absorbed on GN and PyrroN, and is weakly adsorbed on A-pyN-H (-0.20 eV), Z-pyN-H (-0.43 eV), and Z-pyN-O (−0.59 eV) (Table 1). For the remaining structures, the adsorption of COOH is relatively strong, with the adsorption energy ranging from −0.95 to −2.48 eV. However, HCOO exists only on four N-doped graphene structures, that is, A-pyN, A-pyN-O, A-pyN-OH, and Z-GN. The adsorption energies for the four structures are in the range of −1.87 eV∼ −0.84 eV (Table 1).

For HCOOH, the adsorption energies for the studied compounds are in the range of −0.43 ∼ −0.06 eV, which are relatively weak and facilitate its desorption from the catalyst surface. Similar to the HCOOH molecule, the adsorption energies of CO are in the range of −0.26 ∼ −0.03 eV (Table 1).

The possible reaction pathways for the studied compounds are summarized in Figure 2. Based on the computational hydrogen electrode (CHE) model (Norskov et al., 2004), the limiting potential is obtained by U_L = −ΔG_MAX/ e, where ΔG_MAX denotes the maximum free energy difference between the two successive reaction steps. The reduction step corresponding to the limiting potential is defined as the potential determining step (PDS).

As shown in Figures 2A–E, the energy of CO₂ increases by 0.25–0.37 eV from the free molecule to the adsorbed state. After CO₂ is adsorbed on the catalyst surface, it would be hydrogenated by (H⁺ + e⁻) pair. The formation of an O-H bond would produce COOH, while the formation of the C-H bond would generate the HCOO intermediate.

The reaction of CO₂+H⁺+e⁻→*COOH on A-GN, A-pyN, A-pyN-H, A-pyN-O, and A-pyN-OH is uphill by 1.16, 0.84, 2.14, 1.38, and 1.41 eV, respectively. For CO₂+ H⁺+e⁻→*HCOO, the energy increases by 1.73, 1.38, and 1.81 eV for A-pyN, A-pyN-O, and A-pyN-OH, respectively.

The hydrogenation of COOH would generate COHOH, HCOOH, and CO + H₂O. Due to the large energy increase for producing COHOH, that is, 1.55, 0.96, and 1.26 eV for A-GN, A-pyN, and A-pyN-OH, respectively, further discussion is omitted. In COOH, if the OH moiety binds (H⁺+e⁻), it would produce CO + H₂O. If the carbon atom in COOH binds (H⁺+e⁻), it would produce HCOOH. The production of HCOOH and CO is all thermodynamically downhill.

Similarly, the hydrogenation of HCOO may produce H₂COO and HCOOH. As *HCOO→*H₂COO step is endothermic with a large free energy increase (0.88 eV for A-pyN, 0.93 eV for A-pyN-O, and 1.86 eV for A-pyN-OH), further discussion is not provided. Thus, the final product from HCOO is HCOOH.

As illustrated in Figures 2A–E, the COOH intermediate has better performance in producing HCOOH than HCOO. For CO₂→ *CO₂→*COOH→*HCOOH/*CO, the PDS is *CO₂→*COOH (Table 2), which is in agreement with the previous study (Wu et al., 2015). According to the free energy barrier (Figures 2A–E), A-pyN exhibits the highest catalytic activity toward HCOOH with a free energy barrier of 0.84 eV (Figure 2B). The order of catalytic activity for COOH to HCOOH/CO is A-pyN > A-GN > A-pyN-O > A-pyN-OH > A-pyN-H. In addition, CO₂→ *CO₂→*COOH→*CO+*H₂O is the secondary pathway with slightly larger endothermic energy than CO₂→ *CO₂→*COOH→*HCOOH.

The reaction pathways on N-doped zigzag graphene nanoribbons (Figures 2F–I) are similar to those on N-doped armchair graphene nanoribbons. The HCOO intermediate could only stably exist on Z-GN among these N-doped zigzag graphene nanoribbons. To produce HCOOH, the CO₂→*CO₂→*COOH→*HCOOH pathway is more favorable than the CO₂→*CO₂→*HCOO→*HCOOH pathway (Figure 2F). In particular, on Z-GN, the hydrogenation of HCOO generates not only HCOOH but also O + HCHO with an energy barrier of 0.40 eV (Figure 2F). As illustrated in Figure 3, after the formation of O + HCOO, the remaining O atom could be easily hydrogenated into water due to the downhill process. The PDS for producing HCHO is the HCOO formation step with U _L = -0.81 V.

For the *CO₂→*COOH step, it occurred on Z-GN and Z-pyN most easily, in which the energy is uphill by 0.39 eV for Z-GN and downhill by -0.11 eV for Z-pyN, respectively (Figures 2F,G). While for the other two structures, large uphill energy barriers are required, that is, 1.88 eV for Z-pyN-H and 1.72 eV for Z-pyN-O, respectively. After the formation of COOH, its hydrogenation may generate HCOOH, CO + H₂O, or COHOH, in which the formation of HCOOH is the most favorable, followed by CO + H₂O and COHOH. Our calculations indicated that the COOH intermediate on Z-pyN needs an energy barrier of 0.83 eV to form COHOH (Figure 4). After the formation of COHOH, an energy increase of 0.41 eV is required to produce COH + H₂O. The further hydrogenation of COH is relatively easy due to the downhill energy process to release the two competitive final products, that is, CH₃OH and CH₄. A previous study indicated that the formation of CH₄ and CH₃OH is through CO intermediate (Hori et al., 2008), which is different from our results.

As mentioned above, the pyrrolic N-doped structure has no catalytic activity for CO₂RR. For GN, the free energy increase is the largest among all the N-doped graphene structures (2.55 eV). Thus, the catalytic activity of GN is omitted.

In a word, for the studied structures, the most favorable product is HCOOH, followed by CO and COHOH. In particular, the formation of HCOOH and CO is competitive since the free energy of CO is more thermodynamically favorable by only about 0.20 eV than that of HCOOH. This energy difference is similar to the value of 0.28 eV reported earlier (Liu et al., 2016). In a word, Z-pyN and Z-GN possess the highest catalytic activity toward HCOOH due to the smallest limiting potential of −0.38 and −0.39 V, respectively (Table 2), which is lower than −0.44 for PyrroN3 (Liu et al., 2016).

Hydrogen evolution reaction (HER) is the competitive reaction for CO₂RR since the evolution of H would consume the proton–electron pair (H⁺+e⁻) and passivate the catalytic activity of CO₂RR. For the studied structures, the results showed that Z-pyN-O and Z-pyN have large energetic downhill for the adsorption of H*, indicating the enhanced HER in thermodynamic (Figure 5). For Z-GN and A-pyN, they have a negligible free energy barrier (0.03 and 0.04 eV) of H*. For the remaining structures, HER is hindered by large free energy barriers. Therefore, for the most favorable Z-pyN and Z-GN, CO₂RR would be suppressed by HER. However, by choosing a suitable electrolyte, the activation energy of HER would be increased. For instance, according to the expression ΔG _pH= 2.303k _BT pH, in which pH = 0 is selected in the above study, ΔG _pH = 0.42 eV is obtained for pH = 7.0. Thus, the activation energy of HER on Z-GN would be increased from −0.03 to 0.39 eV, comparable to the free energy barrier of 0.38 eV in the CO₂RR process. Thus, the HER could be suppressed by increasing the pH value for Z-GN. While for Z-pyN, CO₂RR is passivated by too strong HER. In a word, Z-GN could be selected as a promising nonmetal electrocatalyst for CO₂RR in generating HCOOH.