Minimum conditions for accurate modeling of urea production via co-electrolysis

Co-electrolysis of carbon oxides and nitrogen oxides promise to simultaneously help restore the balance of the C and N cycles while producing valuable chemicals such as urea. However, co-electrolysis processes are still largely inefficient and numerous knowledge voids persist. Here, we provide a solid thermodynamic basis for modelling urea production via co-electrolysis. First, we determine the energetics of aqueous urea produced under electrochemical conditions based on experimental data, which enables an accurate assessment of equilibrium potentials and overpotentials. Next, we use density functional theory (DFT) calculations to model various co-electrolysis reactions producing urea. The calculated reaction free energies deviate significantly from experimental values for well-known GGA, meta-GGA and hybrid functionals. These deviations stem from errors in the DFT-calculated energies of molecular reactants and products. In particular, the error for urea is approximately -0.25 ± 0.10 eV. Finally, we show that all these errors introduce large inconsistencies in the calculated free-energy diagrams of urea production via co-electrolysis, such that gas-phase corrections are strongly advised.

T he electrochemical co-reduction of species containing nitrogen and carbon to produce chemical commodities can be carried out using renewable electricity [1][2][3][4] , simultaneously aiding to restore the severely imbalanced cycles of nitrogen and carbon [5][6][7] .Among the potential products, urea (CO(NH 2 ) 2 ) is an appealing C-N compound given its enormous relevance in modern agriculture [8][9][10] , and the large energy demands for its industrial production 2,9,11,12 .Although electrocatalytic urea production from N-and C-oxides has been studied at the laboratory scale for more than two decades [13][14][15][16][17][18] , critical challenges regarding the electrochemistry of the C-N coupling are yet to be solved before industrial applications are at hand 19,20 .Some of these challenges are the large associated overpotentials, elusive reaction mechanisms, and low selectivity caused by the concurrent formation of H 2 , CO, formic acid (HCOOH), ammonia (NH 3 ), and other single-carbon and/or single-nitrogen species 16,[21][22][23] .
Either in tandem with experiments or in standalone computational studies, density functional theory (DFT) methods have been extensively used to investigate the electrosynthesis of urea from N-and C-containing feedstocks and design enhanced catalysts 2,[24][25][26][27][28] .These studies frequently use exchangecorrelation (xc) functionals following the generalized gradient approximation (GGA), as they provide a reasonable tradeoff between computational cost and accuracy for the properties of molecules and surfaces [29][30][31][32] .A recent example is the work of Wan et al., in which DFT-based thermodynamic and kinetic models were proposed to explain the selective C-N bond formation on Cu electrodes and rationalize the experimental observations of Shibata et al. 14,16,33 using the BEEF-vdW functional 34 .Notwithstanding, the limitations of GGA functionals are well known for describing gaseous molecules containing multiple bonds, such as O 2 30,35,36 , N 2 and NO x [37][38][39] , and carbon-containing species [40][41][42][43] .These limitations can cause large discrepancies between calculated and experimental equilibrium potentials and impair the predictive capabilities of GGA-based heterogenous (electro)catalytic models, where molecules and surfaces are simultaneously involved 39,44,45 .
Approaches to overcome some of the shortcomings of GGA functionals, such as meta-GGA 46 and hybrid functionals 47 , tend to perform better for the prediction of gas-phase thermochemistry.Unlike GGA functionals, meta-GGA functionals include an approximate dependence on the kinetic energy density 48 , while hybrid functionals incorporate a proportion of exact nonlocal Fock exchange [49][50][51][52] .Interestingly, previous works have shown that when GGA, meta-GGA, and/or hybrid functionals are used to model various families of C- 40,43 and N-containing compounds 38,39,53 , H 2 O 2(g) and O 2(g) 36,44,45 , sizable gas-phase errors are still found.Such errors are systematic and can be mitigated by means of inexpensive semiempirical corrections 38,40,41,43,53 .This strongly suggests that a cautious and early assessment of gas-phase errors is needed to guarantee the accuracy of (electro)catalytic models based upon DFT calculations.
Herein, we study the co-electrolysis of different nitrogen (N 2(g) , NO (g) , NO À 3 aq ð Þ ) and carbon oxides (CO (g) , CO 2(g) ) as feedstocks to produce aqueous urea (CO(NH 2 ) 2(aq) ) using several exchangecorrelation functionals: four GGA functionals, two meta-GGA functionals, and two hybrid functionals.For most gas-phase compounds under study at these three levels of DFT, we pinpoint and correct large gas-phase errors in the calculated energies.Our results stress the importance of gas-phase error assessment in computational electrocatalysis and provide an accurate starting point for modeling urea production on real catalysts by coelectrolysis of CO x and NO x feedstocks.

Methodology
Computational methods.The Vienna ab initio simulation package (VASP) 54 was used to perform the DFT calculations of H 2(g) , N 2(g) , O 2(g) , H 2 O (g) , NH 3(g) , CO(NH 2 ) 2(g) , CO (g) , CO 2(g) , NO (g) , HNO 3(g) , and C (s) .All compounds were modeled in their gas-phase in boxes of 15 × 15 × 15 A ̊3 (in some cases, we changed the size of the vectors by ±0.1 Å to see if more negative energies were found, which was the case only for NO).C (s) was represented here by graphene as a reasonable DFT model of graphite.The latter approximation is enabled by the fact that the interlayer cohesive energy of graphite is small (0.031-0.064 eV) [55][56][57][58][59] .The calculations were carried out for a range of DFT functionals ascending the so-called "Jacob's ladder" 60 : namely GGAs (PBE 61 , PW91 62 , RPBE 63 , BEEF-vdW 64 ), meta-GGAs (TPSS 48 , SCAN 65 ), and hybrids (PBE0 66 , B3LYP 67 ).The C-C distances for graphene obtained in all cases were close to the experimental value of 1.42 A ̊(PBE: The projector augmented-wave (PAW) method was used to represent the interactions between core electron density and valence electrons 69 .A plane-wave cutoff of 450 eV was used in all calculations, assuring converged ΔZPE and reaction energies for the gaseous urea production from N 2(g) and CO 2(g) using PBE In fact, the difference between the reaction energy and ΔZPE obtained with this cutoff differed only by ~0.01 eV from those obtained with a tighter cutoff of 1000 eV (see Supplementary Fig. 1 and Supplementary Table 1 in Supplementary Note 1).The geometry of each molecule was relaxed using the conjugate gradient algorithm until the final forces between the atoms were lower than 0.01 eV A ̊−1 .Gaussian smearing with an electronic temperature of 10 −3 eV was used to ease the convergence of the self-consistent field procedure and, upon convergence, all energies were extrapolated to 0 K. Since the code used is intrinsically periodic, the calculations for molecules were carried out at the Γ-point.Conversely, for graphene a Monkhorst-Pack grid 70 of 8 × 8 × 1 special k-points was used.Spin-unrestricted calculations were performed for O 2(g) (triplet) and NO (g) (doublet).Further details of the input files used to perform the calculations are provided in Supplementary Note 7 and the coordinates of the converged geometries are given in Supplementary Note 8.
The thermodynamic analyses in this study are based upon free energies.The free energy of compound i (G DFT i ) is approximated by means of its DFT energy (E DFT i ), the calculated zero-point energy (ZPE i ) using harmonic frequencies, the difference between the formation enthalpies at 298.15 and 0 K (Δ f H i @298:15K À Δ f H i@0K ), and the entropic contributions (TS i ) taken from thermodynamic tables at T = 298.15K, as shown in Eq. 1 [71][72][73][74][75] .Supplementary Table 2 compiles the DFT-calculated energies, ZPEs, experimental TS, and the differences between the experimental formation enthalpies between 0 and 298.15K for the compounds under study (see also Supplementary Table 5).We provide more details of the thermal contributions in Supplementary Note 4.
Co-electrolysis modeling.We consider the six reactions shown below in which urea is produced by the simultaneous reduction of different C-and N-containing species.
These reactions involve reactants and products in different physical states, such that, under the appropriate external potential, gaseous and aqueous compounds react to produce hydrated urea and liquid water.As DFT simulations of liquids and aqueous systems are possible but challenging and time-consuming and additional statistical analyses are necessary, the DFT values of the corresponding gasphase references, which are rapidly calculated, serve as the basis to estimate their energetics via semiempirical considerations, as depicted in Fig. 1.Moreover, as computational solvation methods have discrepancies with respect to experiments 76 that are larger than the accuracies of the experimental measurements, we do not expect that calculating the solvation energies of the species using DFT would yield lower errors.
The scheme in Fig. 1 shows the energy differences between the states of a generic compound HX. Figure 1 along with the additional considerations detailed below were used to calculate the free energies of all the compounds in the co-electrolysis reactions.We note that thermodynamic cycles based on experimental equilibrium potentials have also been used to semiempirically obtain the free energies of ionic species from the DFT-calculated energies of neutral solids and gaseous compounds 37,77,78 .In the Supplementary Note 6, we show in a stepwise fashion how the free energy of nitrate can be estimated using this approach.
(i) The energetics of proton-electron pairs was calculated by means of the computational hydrogen electrode, which is based on the following equilibrium in solution: 79 .(ii) Based on Fig. 1, the free energy of formation of aqueous urea Þ was estimated by adding the experimental . The experimental solvation energy can be obtained as the difference between the experimental formation energy of aqueous urea (iii) Following Fig. 1, the free energy of formation of liquid water (Δ f G 0 ) was obtained by subtracting the experimental water vaporization energy ðΔ vap G exp H 2 O ¼ 0:09eVÞ 72 from the DFT-calculated free energy of formation of water in the gas phase ðμ DFT 37 , here we use 1 2 H 2(g) and HNO 3(g) as references, as shown in Eq. 8: We note that HNO 3(g) dissociation in Eq. 8 is complete, as it is a strong acid 80 .The free energy of Eq. 8 ðΔG exp 8 Þ can be expressed Fig. 1 Thermodynamic framework.Scheme relating the energy differences between the states of a generic compound HX.The subindices s, l, and g represent HX in the solid, liquid, and gas phases, respectively.The total energy of HX in the gas phase (in green) can be estimated from DFT, using Eq. 1.
The subindex aq refers to hydrated HX, i.e, HX (s) surrounded by water.In red, an anion is produced from the dissociation of HX (aq) .Δ sol G is the free energy of solution (the energy associated to the dissolution of one mole of HX (s) in an infinite amount of water); Δ fus G is the fusion free energy; Δ vap G is the vaporization energy; Δ solv G is the solvation free energy, defined as the energy required to bring a mole of HX (g) from vacuum to a water reservoir; Δ diss G is the dissociation free energy in solution.The dashed red lines indicate that dissociation occurs only for HNO 3 in this study. as: As shown in Fig. 1, ΔG exp 8 corresponds to the sum of the solvation and dissociation energies of HNO 3 ðΔ solvþdiss G exp HNO 3 Þ, which can be calculated as the difference between the experimental formation energy of NO À 3 aq ð Þ (−1.15 eV) and that of HNO 3(g) (−0.76 eV) 72 .Hence, in this study ΔG exp 8 ¼ Δ solvþdiss G exp HNO 3 ¼ À0:39eV.Based on Eq. 9, the free energy of formation of NO À 3 aq ð Þ can be assessed from DFT and experimental data as: where the energetics of protons (μ DFT H þ ) is obtained by invoking the computational hydrogen electrode 79 .With these considerations, DFT and experimental values can be combined to semiempirically calculate the free energies of each compound in Eqs. 2 to 7, and the associated free energies of reaction.
Gas-phase error assessment.An important consideration in the modeling of heterogenous (electro)catalytic reactions is the detection and correction of the errors in the DFT-calculated energies of gas-phase compounds.The errors of O 2(g) , N 2(g) , NO (g) , HNO 3(g) , CO (g) , and CO 2(g) , have previously been calculated using several functionals and large values have been reported in various cases 36,[38][39][40]45,53 . Thesesignificant errors prevent accurate estimations of important quantities in catalysis for the three reasons detailed below.
First, an accurate equilibrium potential for a given reaction may only be rationally obtained by correcting the gas-phase errors of all reactants and products.This is because the equilibrium potential (U eq ) is a function of the reaction free energy (Δ r G) and the number of electrons transferred (n).For a reduction reaction: U eq ¼ ÀΔ r G=n: For example, the reaction in Eq. 2 has an experimental Δ r G of -0.69 eV and involves 4 protonelectron transfers.Thus, its experimental U eq is À0:69 eV À4e À ¼ 0:17 V. Now, the Δ r G using uncorrected PBE is -1.68 eV and U eq is À1:68 eV À4e À ¼ 0:42 V, which deviates by 0.250 V from the experimental value.When the errors of N 2 and CO (the reactants) are corrected, the new PBE equilibrium potential is 0.22 V, which is 0.05 V away from the experimental value.After correcting the error in urea, the PBE calculations match the experimental value.
Second, if the potential-limiting step involves molecules, correcting the gas-phase error or not may lead to different qualitative and quantitative conclusions because the reaction energy experiences a shift.
Third, when the overpotential is calculated, gas-phase errors are always important because the overpotential is the difference between a given potential and the one at equilibrium (for a reduction reaction: η ¼ U eq À U).
It has also been shown that gas-phase errors can affect adsorption-energy scaling relations and volcano plots 36,39,44 , impairing the predictive capability of descriptor-based models of customary use in computational electrocatalysis.
To obtain a general expression to assess the gas-phase errors, we first consider the formation reaction of a hypothetical compound where α; β; γ; and δ are integers, the molecules are in their standard states, and C s ð Þ is modeled as graphene.The total error in the DFT-calculated free energy of formation of as shown in Eq. 12: 36,[38][39][40]44,45,53 ε In addition, the total error is the difference between the individual errors of the products and reactants of Eq. 11: 36,[38][39][40]44,45,53 The DFT error of N 2(g) is calculated from the ammonia synthesis reaction ( 12 35,79 .Assuming that DFT provides an accurate description of the energetics of H 2(g) and C (s 30 , combining Eqs. 12 and 13, and reorganizing, we find an expression for the assessment of the gas-phase error of H α C β N γ O δ , see Eq. 14. A detailed example of the use of these equations for urea is presented in Supplementary Note 2. We note that the experimental standard free energies of formation and heats of formation of the molecules under study are reported within chemical accuracy (4.18 kJ mol −1 or 0.04 eV).In fact, the errors reported on the NIST website for the heats of formation of CO 2 , CO, water, and urea are 0.13, 0.17, 0.040, and 1.2 kJ mol −1 (refs. 71,81), respectively, and the errors are smaller in the ATcT database 82 .In addition, the FreeSolv database commonly reports 2.51 kJ mol −1 as the uncertainty of the experimental hydration free energies 76 .Hence, the experimental values of interest are known to a greater precision compared to the DFT results, which usually involve errors above 0.1 eV (~10 kJ mol −1 ).
We remark that errors may also exist in the adsorbed state.However, as GGA functionals accurately describe the atomic structure, cohesive energy, and bulk moduli of transition metals and their low Miller-index surfaces 29,31 , we expect these errors to be smaller than those of the gas phase.The error assessment detailed in this section may be extended to adsorbates if accurate experimental adsorption energies are available, but these are scarce in the literature 83 .We are aware of two approaches to estimate errors in the adsorbed state.Based on uncertainty considerations, the first method links gas-phase errors to those of the corresponding adsorbates and provides specific corrections for a given species on a substrate (e.g, *COOH on Cu(111)) 43 .Without comparing to experiments, the second method identifies systematic errors of a given adsorbate on a substrate by comparing the adsorption energies using a variety of DFT setups, (i.e, *O on RuO 2 (110)) 84 .
To close this section, we remark that including or neglecting the thermal contributions in Eq. 1 may shift the values of the gasphase errors.These effects are detailed in Supplementary Note 4. Supplementary Table 6 shows the difference of the errors with and without thermal enthalpic contributions from 0 to 298.15 K.

Results
Errors in co-electrolysis reactions.The experimental and DFTcalculated energies of the studied reactions are shown in Supplementary Table 4 and Fig. 2.These values were obtained using the uncorrected DFT energies in Supplementary Table 3 as explained in Supplementary Note 3.For each functional used, the mean absolute errors (xc-MAEs) and maximum absolute errors (xc-MAXs) with respect to experiments are shown in Supplementary Table 4.In addition, Supplementary Table 4 contains the mean and maximum absolute errors for each reaction (r-MAE and r-MAX, respectively).
The r-MAEs and r-MAXs in Supplementary Table 4 indicate that, for the six co-electrolysis reactions studied here, pure GGAbased DFT calculations yield significant deviations with respect to experiments, with r-MAEs spanning from 0.47 to 1.47 eV, and r-MAX values in the range of 1.00 to 2.64 eV.Note in passing that the reaction energies, r-MAE and r-MAX tend to increase as the reactants are more oxidized, as shown in Supplementary Fig. 2.Moreover, the large mean and maximum absolute errors for each functional in Supplementary Table 4 (xc-MAE and xc-MAX), indicate that none of the studied functionals correctly describes all the co-electrolysis reactions, regardless of the functional rung on Jacob's ladder.In fact, for GGA functionals the average xc-MAE and xc-MAX are 0.97 and 2.02 eV; for meta-GGA functionals they are 0.88 and 1.80 eV; and for the hybrid functionals they are 0.49 and 0.84 eV, respectively.Hence, there is an error decrease upon climbing Jacob's ladder, but even hybrid functionals display considerable deviations.
The panels in Fig. 2 aid in visualizing the large discrepancies between theory and experiments.In all the reactions, most of the bars lie far from the experimental value (dashed line in Fig. 2).Consistent with previous works, PBE and PW91 display comparable errors for all reactions 85,86 .Interestingly, BEEF-vdW, RPBE, and TPSS display similar reaction energies in all panels of Fig. 2. In contrast, SCAN and TPSS present large differences although both are meta-GGA functionals.Some similarities are observed in panels e and f for the hybrids, where PBE0 performs better than B3LYP, but significant differences are observed as the N-containing reactant becomes less oxidized Fig. 2 Errors in co-electrolysis reactions.Free energies (Δ r G) of six co-electrolysis reactions calculated with several exchange-correlation functionals: urea production from the co-electrolysis of (a) and CO 2 g ð Þ .In all panels, the respective experimental energy is shown as a dashed black line.The DFT energies of the molecules do not include any gasphase corrections, see the values of these corrections in Table 1.
(panels a-d) and the PBE0 accuracy worsens with respect to B3LYP.As expected, the two hybrid functionals provide more accurate values than GGA and meta-GGA functionals.Overall, calculations using B3LYP lead to results with the smallest errors, presumably as a consequence of its parameterization based on thermochemical data such as atomization energies and ionization potentials.However, for some reactions, the errors of the hybrid functionals surpass the accuracy necessary to allow for accurate predictions (<0:1eV): B3LYP yields errors of −0.45, −0.31, and 0.36 eV for reactions a, b, and f, and PBE0 yields errors of −1.23, −1.00, −1.04, −0.81, and −0.28 eV for the reactions in Fig. 2a-e.
Errors in the molecules.Table 1 summarizes the gas-phase errors of the species involved in the co-electrolysis reactions for all the functionals under analysis using Eq.14.Note that NO À 3 aq ð Þ and CO(NH 2 ) 2(aq) display the same errors as their respective gaseous references, HNO 3(g) and CO(NH 2 ) 2(g) .This is because the energies of these species were calculated semiempirically from DFT energies of the gases and experimental values (see section 2.2).
As shown in Table 1, the DFT errors of most species under study are significant regardless of the functional rung on Jacob's ladder, and reach in some cases values more negative than −1 eV.This is the case of HNO 3(g) using RPBE (−1.04 eV), BEEF-vdW (−1.26 eV), and TPSS (−1.19 eV).The hybrid functionals B3LYP and PBE0 yield the lowest gas-phase errors (MAEs of 0.17 and 0.24 eV in Table 1) but still exhibit large MAX figures (0.28 and 0.71 eV).In fact, PBE0 displays the largest error for N 2(g) (0.71 eV), which partly explains the substantial deviations of this functional in Fig. 2a, b.
The values in Table 1 can be employed to rapidly estimate the error cancellation of a functional when modeling a chemical reaction.Error cancellation may lead to accurate predictions of reaction energies.For instance, in Fig. 2f we observe for PBE0 an almost complete error cancellation because the errors of reactants and products differ only by 0.05 eV: Þ%À0:05eV.Furthermore, the errors in N 2(g) and O 2(g) are large for all functionals, spanning from -0.16 to 0.71 eV for N 2(g), and from -0.80 to -0.17 eV for O 2(g) .For urea, significant errors are found for all scrutinized functionals with PBE0 presenting the lowest value (-0.11 eV) and, surprisingly, BEEF-vdW displaying the largest (-0.40 eV).It is worth noting that this range of errors for urea is narrow compared to the other molecules in Table 1 and that all values are negative.Hence, a tentative estimate for the DFT-based error for the formation energy of urea is -0.25 ± 0.10 eV, which corresponds to the average and standard deviation of the corresponding values in Table 1.
Figure 3 provides a graphical representation of the values in Table 1.The ranges of the errors are larger than 0.20 eV in all cases.For most molecules, we observe that RPBE, BEEF-vdW, and TPSS exhibit the largest negative errors.In contrast, PBE0 always yields the largest positive deviations for N 2(g) and CO (g) but also the lowest errors for O 2(g) , urea and HNO 3(g) /NO - 3(aq) .Finally, while extent of the error fluctuation of PBE, PW91, and SCAN over the whole set of molecules is rather similar, B3LYP shows the smallest and relatively stable set of errors, all relatively close to zero.However, as mentioned before, our general conclusion is that none of the functionals in Table 1 yields satisfactory energies for the co-electrolysis reactions under analysis.
Importantly, if experimental results are not available to calculate the gas-phase errors using Eq.14, one could rely on highly accurate quantum chemical methods based on wave-function theory such as CCSD(T) using large basis sets.Alternatively, one can use correction approaches based on structural features, such as the number of oxygen atoms in the molecule 39 , the presence of certain functional groups 38 , or the occurrence of specific chemical structures within the compound, such as CO-, OCO-, ONO-, NNO-, or -NOH backbones 41,43,53 .For instance, from a functional group perspective, CO(NH 2 ) 2(g) can be considered an amide with an amino group bound to it.The respective PBE errors of the amino and amide groups are 0.00 and −0.17 eV 38 , yielding a total error of −0.17 eV, which agrees well with that in Table 1 (−0.20 eV).This approximation is somewhat satisfactory for the other GGAs studied (PW91: −0.15 eV, RPBE: −0.16 eV, BEEF-vdW: −0.47 eV 38 vs −0.20, −0.22 and −0.40 eV in this study).However, we note that ð Þ yields a more accurate approximation (PBE: −0.17 eV, PW91: −0.15 eV, RPBE: −0.18 eV, BEEF-vdW: −0.42 eV) 38 .This is because of the double counting of one of the C-N bonds: the amide correction was designed to account for the error in the bond between an sp 3 C and -CO(NH 2 ), and the amine correction for the bond between an sp 3 C and -NH 2 .
Implications for electrocatalysis.To illustrate the effect of gasphase errors on electrocatalysis, the DFT-uncorrected values in Supplementary Table 4 and the errors in Table 1 were used to build free-energy diagrams of the thermodynamically ideal catalyst for each co-electrolysis reaction, see Fig. 4 and Supplementary Figures 3-7 in the Supplementary Note 5.The concept of an ideal catalyst is commonly employed in electrocatalysis to outline the most efficient conversion that conforms to the first and second laws of thermodynamics [87][88][89] , thus serving as a benchmark for real catalysts.In an ideal catalyst, the reaction energy of all electrochemical steps is the same and corresponds to the overall reaction energy divided by the number of electrons transferred (i.e, ΔG i ¼ Δ r G=n).Numerically, the magnitude of the ideal electrochemical steps is identical to the equilibrium potential.Hence, the chemical identity of the intermediates need not be known to build the ideal free-energy diagram.In contrast, real catalysts usually display asymmetric free-energy diagrams and require knowledge of the chemical identity and energetics of the intermediates.For the co-reduction of NO À 3 aq ð Þ and CO 2(g) in Fig. 4, previous works proposed the following mechanism: 20,90,91 *NO 3 reduction to *NO 2 , then coupling with *CO 2 to form *CO 2 NO 2 .Subsequent protonation of *CO 2 NO 2 yields *CO 2 NH 2 , which in turn reduces to *COOHNH 2 in the potential-determining step 92 .*COOHNH 2 reduces to *CONH 2 , which couples to *NO 2 , producing *CONO 2 NH 2 .Finally, *CONO 2 NH 2 is hydrogenated twice to give urea.
Figure 4 shows three energy diagrams of the NO À 3 aq ð Þ and CO 2(g) co-electrolysis reaction to urea on the ideal catalyst (Eq.7).Panel a was built with the uncorrected DFT values, i.e, with no gas-phase corrections.In panel b, the error of CO 2(g) (the C-containing reactant) was accounted for, leaving the DFT energies of both NO À 3 aq ð Þ and CO(NH 2 ) 2(aq) uncorrected.In panel c, the DFT-energies of CO 2(g) and NO À 3 aq ð Þ were corrected, while the energy of CO(NH 2 ) 2(aq) remained uncorrected, except for the black line, in which DFT and experiments coincide.
Figure 4a shows that all functionals diverge from the freeenergy profile of the ideal catalysts (calculated on the basis of experimental values) as more electrochemical steps are considered, reaching a maximum deviation at the last step of the catalytic pathway.This maximum deviation corresponds to the difference between the DFT reaction energy and its experimental counterpart.The telescopic effect in Fig. 4a also occurs in panels b and c, but with nuances introduced by the partial corrections.We note the total error can also be obtained by assessing the difference between the errors of reactants and products.For example, based on the values in Table 1, for Eq. 7 and BEEF-vdW the resulting error is 2 4c the departures of the predicted values from calculations using the various functionals with respect to those from experiments stem from the error in CO(NH 2 ) 2(aq) , as it is the only remaining uncorrected species.In other words, the difference between DFT-based and experimental values for the last reaction step of Fig. 4c is . Moreover, we note that after correcting the error of urea using the respective values in Table 1, all the gas-phase errors are accounted for and the "DFT + corrections" diagram becomes that of the ideal catalyst, which is shown in black in all three panels of Fig. 4.
Analogous diagrams for the other reactions under study are given in Supplementary Figs.3-7.We emphasize that the conclusions drawn from Fig. 4 also hold for these Supplementary figures.Since ε COðNH 2 Þ 2 aq ð Þ < 0 for all the functionals assessed (see Table 1), the DFT-calculated lines are always below the experimental values in the bottom panels of Fig. 4 and Supplementary Figs.3-7.

Conclusions
Simultaneous electrocatalytic reduction of nitrogen and carbon pollutants to produce urea is an appealing alternative to help remediate the colossal imbalances of the nitrogen and carbon cycles.Herein, we showed how experimental data can be coupled with DFT-calculated gas-phase energies to model six coreduction reactions of carbon and nitrogen oxides to produce hydrated urea.The average MAE/MAX values versus experiments are 0.97 eV/2.02eV for GGA functionals (PBE, PW91, RPBE, and BEEF-vdW), while those of meta-GGAs (TPSS and SCAN) are 0.88 eV/1.80 eV, and those of the hybrids (PBE0, B3LYP) are 0.49 eV/0.84 eV.Hence, the use of DFT to model these reactions entails large errors, even for hybrid functionals, indicating that accurate predictions are only attained once the DFT errors of all the molecules under study are corrected.
Moreover, the DFT error in the formation energy of urea spans a relatively narrow range of values, such that -0.25 ± 0.10 eV is a reasonable error estimate for DFT calculations, although the use of specific corrections is always more advisable than an average.The effect of these numerical deviations in catalysis was illustrated for the free-energy diagrams of the ideal electrocatalyst extracted from experimental data for various co-electrolysis reactions.The departures of DFT predictions from the experimental trends are substantial for all functionals.However, we showed that the errors can easily be corrected in a semiempirical manner.
All this hints toward the need for an assessment of gas-phase errors at the early stages of computational electrocatalysis research to avoid potentially inaccurate and misleading conclu- sions, regardless of the chosen rung on Jacob's ladder of density functional approximations.

Fig. 1
07eVÞ73 and the experimental formation energy of solid urea

Fig. 3
Fig. 3 Individual errors of the species.DFT errors (ε i ) for the compounds involved in the co-electrolysis reactions.Circles (•) are for GGA errors, triangles (▼) are for meta-GGA errors, and crosses (×) are for hybrid errors.The vertical bars correspond to the ranges spanned by the functionals.All values are in eV.

Fig. 4
Fig.4Free energy diagrams for the co-electrolysis of NO À 3ðaqÞ and CO2(g) to urea on the ideal catalyst.The diagrams were built using (a) uncorrected DFT energies (Supplementary Table4), (b) the uncorrected DFT energies of NO À

4
Fig.4Free energy diagrams for the co-electrolysis of NO À 3ðaqÞ and CO2(g) to urea on the ideal catalyst.The diagrams were built using (a) uncorrected DFT energies (Supplementary Table4), (b) the uncorrected DFT energies of NO À Fig. 4 Free energy diagrams for the co-electrolysis of NO À 3ðaqÞ and CO2(g) to urea on the ideal catalyst.The diagrams were built using (a) uncorrected DFT energies (Supplementary Table 4), (b) the uncorrected DFT energies of NO À 3 aq ð Þ and CO(NH 2 ) 2(aq) , and the corrected energy of CO 2(g) , (c) the corrected energies of NO À 3 aq ð Þ and CO 2(g) , and the uncorrected DFT energy of CO(NH 2 ) 2(aq) .Ideal values from experiments are in black.The y-axis scale is divided into multiples of the experimental equilibrium potential.

Table 1
Individual gas-phase errors.