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

ARTICLE Minimum conditions for accurate modeling of urea production via co-electrolysis Ricardo Urrego-Ortiz 1,2, Santiago Builes3, Francesc Illas 1, Stefan T. Bromley1,4, Marta Costa Figueiredo 5& Federico Calle-Vallejo 2,6✉ 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 coelectrolysis. 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, metaGGA 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. https://doi.org/10.1038/s42004-023-00990-7 OPEN 1Departament de Ciència de Materials i Química Física & Institut de Química Teòrica i Computacional (IQTCUB), Universitat de Barcelona, C/ Martí i Franquès 1, 08028 Barcelona, Spain. 2Nano-Bio Spectroscopy Group and European Theoretical Spectroscopy Facility (ETSF), Department of Polymers and Advanced Materials: Physics, Chemistry and Technology, University of the Basque Country UPV/EHU, Av. Tolosa 72, 20018 San Sebastián, Spain. 3Escuela de Ciencias Aplicadas e Ingeniería, Universidad EAFIT, Carrera 49 # 7 sur 50, 050022 Medellín, Colombia. 4Institució Catalana de Recerca i Estudis Avançats (ICREA), Passeig Lluís Companys 23, 08010 Barcelona, Spain. 5Eindhoven Institute of Renewable Energy Systems (EIRES), Eindhoven University of Technology, PO Box 513 Eindhoven 5600 MB, The Netherlands. 6IKERBASQUE, Basque Foundation for Science, Plaza de Euskadi 5, 48009 Bilbao, Spain. ✉email: [email protected] COMMUNICATIONS CHEMISTRY | (2023) 6:196 | https://doi.org/10.1038/s42004-023-00990-7 | www.nature.com/commschem 1 1234567890():,; The electrochemical co-reduction of species containing nitrogen and carbon to produce chemical commodities can be carried out using renewable electricity1–4, simultaneously aiding to restore the severely imbalanced cycles of nitrogen and carbon5–7. Among the potential products, urea (CO(NH 2 ) 2 ) is an appealing C-N compound given its enormous relevance in modern agriculture8–10, and the large energy demands for its industrial production2,9,11,12. Although electrocatalytic urea production from Nand C-oxides has been studied at the laboratory scale for more than two decades13–18, critical challenges regarding the electrochemistry of the C-N coupling are yet to be solved before industrial applications are at hand19,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 species16,21–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 Nand C-containing feedstocks and design enhanced catalysts2,24–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 surfaces29–32. A recent example is the work of Wanetal.,inwhichDFT-basedthermodynamicandkinetic 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 functional34. 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–39, and carbon-containing species40–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 involved39,44,45. Approaches to overcome some of the shortcomings of GGA functionals, such as meta-GGA46 and hybrid functionals47, 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 density48, while hybrid functionals incorporate a proportion of exact nonlocal Fock exchange49–52. Interestingly, previous works have shown that when GGA, meta-GGA, and/or hybrid functionals are used to model various families of C40,43 and N-containing compounds38,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 corrections38,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 3aq ðÞ ) 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–59. The calculations were carried out for a range of DFT functionals ascending the so-called “Jacob’s ladder”60: namely GGAs (PBE61, PW9162, RPBE63, BEEF-vdW64), meta-GGAs (TPSS48, SCAN65), and hybrids (PBE066, B3LYP67). The C-C distances for graphene obtained in all cases were close to the experimental value of 1.42 A(PBE: 1.43 A, PW91: 1.43 A, RPBE: 1.43 A, BEEF-vdW: 1.43 A, TPSS: 1.42 A, SCAN: 1.42 A, B3LYP: 1.42 A, PBE0: 1.42 A)68. The projector augmented-wave (PAW) method was used to represent the interactions between core electron density and valence electrons69. 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 (N2g ðÞþCO2g ðÞþ3H2g ðÞ!CO NH2  2g ðÞþH2Og ðÞ ). 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−3eV 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 grid70 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(GDFT i) is approximated by means of its DFT energy (EDFT i), the calculated zero-point energy (ZPEi) using harmonic frequencies, the difference between the formation enthalpies at 298.15 and 0 K (ΔfHi@298:15KΔfHi@0K), and the entropic contributions (TSi) taken from thermodynamic tables at T=298.15 K, as shown in Eq. 171–75. Supplementary Table 2 compiles the DFT-calculated energies, ZPEs, experimental TS, and the differences between the experimental formation enthalpies between 0 and 298.15 K for the compounds under study (see also Supplementary Table 5). We provide more details of the thermal contributions in Supplementary Note 4. GDFT iEDFT iþZPEiþðΔfHi@298:15KΔfHi@0KÞTSið1Þ Co-electrolysis modeling. We consider the six reactions shown below in which urea is produced by the simultaneous reduction of different Cand N-containing species. N2g ðÞþCO g ðÞþ4Hþþ4e!CO NH2  2aq ðÞ ð2Þ N2g ðÞþCO2g ðÞþ6Hþþ6e!CO NH2  2aq ðÞ þH2OlðÞ ð3Þ ARTICLE COMMUNICATIONS CHEMISTRY | https://doi.org/10.1038/s42004-023-00990-7 2COMMUNICATIONS CHEMISTRY | (2023) 6:196 | https://doi.org/10.1038/s42004-023-00990-7 | www.nature.com/commschem 2NO g ðÞþCO g ðÞþ8Hþþ8e!CO NH2  2aq ðÞ þ2H2OlðÞ ð4Þ 2NO g ðÞþCO2g ðÞþ10Hþþ10e!CO NH2  2aq ðÞ þ3H2OlðÞ ð5Þ 2NO 3aq ðÞ þCO g ðÞþ16Hþþ14e!CO NH2  2aq ðÞ þ6H2OlðÞ ð6Þ 2NO 3aq ðÞ þCO2g ðÞþ18Hþþ16e!CO NH2  2aq ðÞ þ7H2OlðÞ ð7Þ 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 experiments76 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. 1shows the energy differences between the states of a generic compound HX. Figure 1along 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 compounds37,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: Hþþe$1 2H2g ðÞ , such that 1 2μ0 H2g ðÞ ¼μ0 Hþþe ðÞ 79. (ii) Based on Fig. 1, the free energy of formation of aqueous urea ðΔfG0 CO NH2 ðÞ 2aq ðÞ Þwas estimated by adding the experimental solvation energy ðΔsolvGexp CO NH2 ðÞ 2aq ðÞ Þto the DFT-calculated formation energy of gaseous urea ðΔfGDFT CO NH2 ðÞ 2g ðÞ Þ, i.e, ΔfG0 CO NH2 ðÞ 2aq ðÞ ¼ΔfGDFT CO NH2 ðÞ 2g ðÞ þΔsolvGexp CO NH2 ðÞ 2aq ðÞ . The experimental solvation energy can be obtained as the difference between the experimental formation energy of aqueous urea ðΔfGexp CO NH2 ðÞ 2aq ðÞ Þand the experimental formation energy of gaseous urea ðΔfGexp CO NH2 ðÞ 2g ðÞ ¼1:57eVÞ71,73.ΔfGexp CO NH2 ðÞ 2aq ðÞ is calculated, in accordance with Fig. 1, by combining the experimental free energy of solution ðΔsolGexp CO NH2 ðÞ 2 ¼0:07eVÞ73 and the experimental formation energy of solid urea ðΔfGexp CO NH2 ðÞ 2sðÞ ¼2:04 eVÞ71,thus ΔfGexp CO NH2 ðÞ 2aq ðÞ ¼2:11eV. Finally, ΔsolvGexp CO NH2 ðÞ 2aq ðÞ ¼ 0:54eV and ΔfG0 CO NH2 ðÞ 2aq ðÞ ¼ΔfGDFT CO NH2 ðÞ 2g ðÞ 0:54eV. (iii) Following Fig. 1, the free energy of formation of liquid water (ΔfG0 H2OðlÞ ) was obtained by subtracting the experimental water vaporization energy ðΔvapGexp H2O¼0:09eVÞ72 from the DFT-calculated free energy of formation of water in the gas phase ðμDFT H2Og ðÞ Þ, i.e, ΔfG0 H2OlðÞ ¼ΔfGDFT H2OðgÞ 0:09eV. (iv) Because calculating the energies of dissolved nitrate ðNO 3aq ðÞ Þ with DFT is problematic37,hereweuse1 2H 2(g) and HNO 3(g) as references, as shown in Eq. 8: HNO3ðgÞ!NO 3ðaqÞþHþð8Þ We note that HNO 3(g) dissociation in Eq. 8is complete, as it is a strong acid80. The free energy of Eq. 8ðΔGexp 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) .ΔsolGis the free energy of solution (the energy associated to the dissolution of one mole of HX (s) in an infinite amount of water); ΔfusGis the fusion free energy; ΔvapGis the vaporization energy; ΔsolvGis the solvation free energy, defined as the energy required to bring a mole of HX (g) from vacuum to a water reservoir; ΔdissGis the dissociation free energy in solution. The dashed red lines indicate that dissociation occurs only for HNO 3 in this study. COMMUNICATIONS CHEMISTRY | https://doi.org/10.1038/s42004-023-00990-7 ARTICLE COMMUNICATIONS CHEMISTRY | (2023) 6:196 | https://doi.org/10.1038/s42004-023-00990-7 | www.nature.com/commschem 3 as: ΔGexp 8¼ΔfG0 NO 3ðaqÞþ1 2ΔfG0 H2ðgÞΔfG0 HNO3ðgÞð9Þ As shown in Fig. 1,ΔGexp 8corresponds to the sum of the solvation and dissociation energies of HNO 3 ðΔsolvþdissGexp HNO3Þ, which can be calculated as the difference between the experimental formation energy of NO 3aq ðÞ(−1.15 eV) and that of HNO 3(g) (−0.76 eV)72. Hence, in this study ΔGexp 8¼ ΔsolvþdissGexp HNO3¼0:39eV. Based on Eq. 9, the free energy of formation of NO 3aq ðÞ can be assessed from DFT and experimental data as: ΔfGDFT NO 3aq ðÞ ¼ΔfGDFT HNO3g ðÞ μDFT Hþ0:39 eV ð10Þ where the energetics of protons (μDFT Hþ) is obtained by invoking the computational hydrogen electrode79. With these considerations, DFT and experimental values can be combined to semiempirically calculate the free energies of each compound in Eqs. 2to 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 cases36,38–40,45,53. These significant 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 (ΔrG) and the number of electrons transferred (n). For a reduction reaction: Ueq ¼ΔrG=n:For example, the reaction in Eq. 2has an experimental ΔrGof -0.69 eV and involves 4 protonelectron transfers. Thus, its experimental Ueq is 0:69 eV 4e¼0:17 V. Now, the ΔrGusing uncorrected PBE is -1.68 eV and Ueq is 1:68eV 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: η¼Ueq U). It has also been shown that gas-phase errors can affect adsorption-energy scaling relations and volcano plots36,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 HαCβNγOδ: α 2H2g ðÞþβCsðÞþγ 2N2g ðÞþδ 2O2g ðÞ!HαCβNγOδð11Þ where α;β;γ;and δare integers, the molecules are in their standard states, and CsðÞ is modeled as graphene. The total error in the DFT-calculated free energy of formation of HαCβNγOδ, denoted εT, is determined as the difference between the DFT prediction (ΔfGDFT HαCβNγOδ) and the experimental value ðΔfGexp HαCβNγOδÞas shown in Eq. 12:36,38–40,44,45,53 εT¼ΔfGDFT HαCβNγOδΔfGexp HαCβNγOδð12Þ In addition, the total error is the difference between the individual errors of the products and reactants of Eq. 11:36,38–40,44,45,53 εT¼εHαCβNγOδα 2εH2g ðÞ βεCsðÞ γ 2εN2g ðÞ δ 2εO2g ðÞ ð13Þ The DFT error of N 2(g) is calculated from the ammonia synthesis reaction (1 2N2g ðÞþ3 2H2g ðÞ!NH3g ðÞ )38 and εO2g ðÞ from the water formation reaction (H2g ðÞþ1 2O2g ðÞ!H2Og ðÞ )35,79. Assuming that DFT provides an accurate description of the energetics of H 2(g) and C (s) (i.e, εH2g ðÞ εCsðÞ 0)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. εHαCβNγOδg ðÞ ¼ΔfGDFT HαCβNγOδg ðÞþγ 2εN2þδ 2εO2  ΔfGexp HαCβNγOδ ð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−1or 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 database82. In addition, the FreeSolv database commonly reports 2.51 kJ mol−1as the uncertainty of the experimental hydration free energies76. 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 surfaces29,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 literature83. 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. 1may 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 ARTICLE COMMUNICATIONS CHEMISTRY | https://doi.org/10.1038/s42004-023-00990-7 4COMMUNICATIONS CHEMISTRY | (2023) 6:196 | https://doi.org/10.1038/s42004-023-00990-7 | www.nature.com/commschem 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 xcMAE 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. 2aid 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 reactions85,86. Interestingly, BEEFvdW, 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 (ΔrG) of six co-electrolysis reactions calculated with several exchange-correlation functionals: urea production from the co-electrolysis of (a)N2g ðÞand CO g ðÞ ,(b)N2g ðÞand CO2g ðÞ ,(c)NO g ðÞand CO g ðÞ ,(d)NO g ðÞand CO2g ðÞ ,(e)NO 3ðaqÞand +CO g ðÞ ,(f)NO 3ðaqÞ and CO2gðÞ . 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. COMMUNICATIONS CHEMISTRY | https://doi.org/10.1038/s42004-023-00990-7 ARTICLE COMMUNICATIONS CHEMISTRY | (2023) 6:196 | https://doi.org/10.1038/s42004-023-00990-7 | www.nature.com/commschem 5 (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 1summarizes 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 3aq ðÞ 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). AsshowninTable1, the DFT errors of most species under study are significant regardless of the functional rung on Jacob’sladder, 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 1can 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: εurea 2εNO 3aq ðÞ εCO2¼ 0:11eV 20:12eVðÞ0:18eVðÞ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 1and 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 3provides 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) /NO3(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 1yields 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 molecule39, the presence of certain functional groups38, or the occurrence of specific chemical structures within the compound, such as CO-, OCO-, ONO-, NNO-, or -NOH backbones41,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 eV38, 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, BEEFvdW: −0.47 eV38 vs −0.20, −0.22 and −0.40 eV in this study). However, we note that εCOðNH2Þ2g ðÞ εamide g ðÞþ1 2εamine g ðÞ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 sp3Cand -CO(NH 2 ), and the amine correction for the bond between an sp3 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 1were used to build free-energy diagrams of the thermodynamically ideal catalyst for each co-electrolysis reaction, see Fig. 4and 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 thermodynamics87–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 Table 1 Individual gas-phase errors. Species PBE PW91 RPBE BEEF-vdW TPSS SCAN B3LYP PBE0 N 2(g) 0.49 0.52 0.11 −0.16 0.00 0.47 0.28 0.71 O 2(g) −0.42 −0.27 −0.70 −0.78 −0.80 −0.40 −0.28 −0.17 NO (g) 0.04 0.15 −0.29 −0.47 −0.41 0.05 0.00 0.26 CO (g) 0.30 0.30 −0.04 −0.13 −0.05 0.30 −0.03 0.40 CO 2(g) −0.14 −0.11 −0.42 −0.52 −0.47 −0.03 −0.18 0.18 HNO 3(g) −0.88 −0.79 −1.04 −1.26 −1.19 −0.51 −0.19 −0.12 NO 3aqðÞ −0.88 −0.79 −1.04 −1.26 −1.19 −0.51 −0.19 −0.12 CO(NH 2 ) 2(g) −0.20 −0.20 −0.22 −0.40 −0.36 −0.32 −0.20 −0.11 CO(NH 2 ) 2(aq) −0.20 −0.20 −0.22 −0.40 −0.36 −0.32 −0.20 −0.11 MAE 0.39 0.37 0.45 0.60 0.54 0.32 0.17 0.24 MAX 0.88 0.79 1.04 1.26 1.19 0.52 0.28 0.71 DFT errors in the formation energy of the species involved in the co-electrolysis reactions. The MAEs and MAXs are reported for each functional. All values are in eV. ARTICLE COMMUNICATIONS CHEMISTRY | https://doi.org/10.1038/s42004-023-00990-7 6COMMUNICATIONS CHEMISTRY | (2023) 6:196 | https://doi.org/10.1038/s42004-023-00990-7 | www.nature.com/commschem reaction energy divided by the number of electrons transferred (i.e, ΔGi¼ΔrG=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 3aq ðÞ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 step92.*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 4shows three energy diagrams of the NO 3aq ðÞ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 3aq ðÞ and CO(NH 2 ) 2(aq) uncorrected. In panel c, the DFT-energies of CO 2(g) and NO 3aq ðÞ 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. 7and BEEF-vdW the resulting error is 2 εNO 3aq ðÞ þεCO2ðgÞεCOðNH2Þ2aq ðÞ ¼ 21:26 0:52 þ0:40 ¼2:64eV. In Fig. 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. 4cis εCOðNH2Þ2aq ðÞ . 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. 4also hold for these Supplementary figures. Since εCOðNH2Þ2aq ðÞ < 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. 4and 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.02 eV 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 concluFig. 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. COMMUNICATIONS CHEMISTRY | https://doi.org/10.1038/s42004-023-00990-7 ARTICLE COMMUNICATIONS CHEMISTRY | (2023) 6:196 | https://doi.org/10.1038/s42004-023-00990-7 | www.nature.com/commschem 7 sions, regardless of the chosen rung on Jacob’s ladder of density functional approximations. Data availability The authors declare that data supporting the findings of this study are available within the paper and its supplemental material file. Additional data are available from the corresponding author upon reasonable request. Received: 27 March 2023; Accepted: 22 August 2023; References 1. Daiyan, R., MacGill, I. & Amal, R. Opportunities and challenges for renewable power-to-X. ACS Energy Lett. 5, 3843–3847 (2020). 2. Chen, C. et al. Coupling N2 and CO2 in H2O to synthesize urea under ambient conditions. Nat. Chem. 12, 717–724 (2020). 3. Jouny, M. et al. Formation of carbon–nitrogen bonds in carbon monoxide electrolysis. Nat. Chem. 11, 846–851 (2019). 4. Chen, L. & Shi, J. Co-electrolysis toward value-added chemicals. Sci. China Mater. 65,1–9 (2022). 5. Zeng, Y., Priest, C., Wang, G. & Wu, G. Restoring the nitrogen cycle by electrochemical reduction of nitrate: progress and prospects. Small Methods 4, 2000672 (2020). 6. Gruber, N. & Galloway, J. N. An earth-system perspective of the global nitrogen cycle. Nature 451, 293–296 (2008). 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 3aqðÞ and CO(NH 2 ) 2(aq) , and the corrected energy of CO 2(g) ,(c) the corrected energies of NO 3aqðÞ 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. ARTICLE COMMUNICATIONS CHEMISTRY | https://doi.org/10.1038/s42004-023-00990-7 8COMMUNICATIONS CHEMISTRY | (2023) 6:196 | https://doi.org/10.1038/s42004-023-00990-7 | www.nature.com/commschem 7. Rockström, J. et al. A safe operating space for humanity. Nature 461, 472–475 (2009). 8. Barzagli, F., Mani, F. & Peruzzini, M. From greenhouse gas to feedstock: formation of ammonium carbamate from CO2 and NH3 in organic solvents and its catalytic conversion into urea under mild conditions. Green Chem. 13, 1267 (2011). 9. Giddey, S., Badwal, S. P. S. & Kulkarni, A. Review of electrochemical ammonia production technologies and materials. Int. J. Hydrog. Energy 38, 14576–14594 (2013). 10. Wang, M. et al. Can sustainable ammonia synthesis pathways compete with fossilfuel based haber–bosch processes. Energy Environ. Sci. 14, 2535–2548 (2021). 11. Glibert, P. M., Harrison, J., Heil, C. & Seitzinger, S. Escalating worldwide use of urea –a global change contributing to coastal eutrophication. Biogeochemistry 77, 441–463 (2006). 12. Meessen, J. Urea synthesis. Chem. Ing. Tech. 86, 2180–2189 (2014). 13. Shibata, M., Yoshida, K. & Furuya, N. Electrochemical synthesis of urea at gasdiffusion electrodes i. simultaneous reduction of carbon dioxide and nitrite ions at Zn catalysts. Denki Kagaku Oyobi Kogyo Butsuri Kagaku 64, 1068–1073 (1996). 14. Shibata, M., Yoshida, K. & Furuya, N. Electrochemical synthesis of urea at gasdiffusion electrodes: ii. simultaneous reduction of carbon dioxide and nitrite ions at Cu, Ag and Au catalysts. J. Electroanal. Chem. 442,67–72 (1998). 15. Shibata, M., Yoshida, K. & Furuya, N. Electrochemical synthesis of urea at gas‐ diffusion electrodes: III. Simultaneous reduction of carbon dioxide and nitrite ions with various metal catalysts. J. Electrochem. Soc. 145, 595–600 (1998). 16. Shibata, M., Yoshida, K. & Furuya, N. Electrochemical synthesis of urea at gas‐ diffusion electrodes: iv. simultaneous reduction of carbon dioxide and nitrate ions with various metal catalysts. J. Electrochem. Soc. 145, 2348–2353 (1998). 17. Shibata, M., Yoshida, K. & Furuya, N. Electrochemical synthesis of urea at gasdiffusion electrodes v. simultaneous reduction of carbon dioxide and nitrite ions with various boride catalysts. Denki Kagaku Oyobi Kogyo Butsuri Kagaku 66, 584–589 (1998). 18. Shibata, M. & Furuya, N. Electrochemical synthesis of urea at gas-diffusion electrodes: vi. simultaneous reduction of carbon dioxide and nitrite ions with various metallophthalocyanine catalysts. J. Electroanal. Chem. 507,177–184 (2001). 19. Li, J., Zhang, Y., Kuruvinashetti, K. & Kornienko, N. Construction of C–N bonds from small-molecule precursors through heterogeneous electrocatalysis. Nat. Rev. Chem. 6, 303–319 (2022). 20. Jiang, M. et al. Review on electrocatalytic coreduction of carbon dioxide and nitrogenous species for urea synthesis. ACS Nano 17, 3209–3224 (2023). 21. Chen, C., He, N. & Wang, S. Electrocatalytic C–N coupling for urea synthesis. Small Sci. 1, 2100070 (2021). 22. Kim, J. E., Choi, S., Balamurugan, M., Jang, J. H. & Nam, K. T. Electrochemical C–N bond formation for sustainable amine synthesis. Trends Chem. 2, 1004–1019 (2020). 23. Mei, Z. et al. Recent progress in electrocatalytic urea synthesis under ambient conditions. ACS Sustain. Chem. Eng. 10, 12477–12496 (2022). 24. Jiao, D. et al. Mo 2 P monolayer as a superior electrocatalyst for urea synthesis from nitrogen and carbon dioxide fixation: a computational study. Energy Environ. Mater.https://doi.org/10.1002/eem2.12496 (2022). 25. Zhu, X., Zhou, X., Jing, Y. & Li, Y. Electrochemical synthesis of urea on MBenes. Nat. Commun. 12, 4080 (2021). 26. Yang, G.-L. et al. Gaseous CO 2 coupling with N-containing intermediates for key C–N bond formation during urea production from coelectrolysis over Cu. ACS Catal. 12, 11494–11504 (2022). 27. Huang, Y. et al. Direct electrosynthesis of urea from carbon dioxide and nitric oxide. ACS Energy Lett. 7, 284–291 (2022). 28. Liu, X., Jiao, Y., Zheng, Y., Jaroniec, M. & Qiao, S.-Z. Mechanism of C-N bonds formation in electrocatalytic urea production revealed by ab initio molecular dynamics simulation. Nat. Commun. 13, 5471 (2022). 29. Janthon, P. et al. Bulk properties of transition metals: a challenge for the design of universal density functionals. J. Chem. Theory Comput. 10, 3832–3839 (2014). 30. Kurth, S., Perdew, J. P. & Blaha, P. Molecular and solid-state tests of density functional approximations: LSD, GGAs, and Meta-GGAs. Int. J. Quant. Chem. 75, 889–909 (1999). 31. Vega, L., Ruvireta, J., Viñes, F. & Illas, F. Jacob’s ladder as sketched by escher: assessing the performance of broadly used density functionals on transition metal surface properties. J. Chem. Theory Comput. 14, 395–403 (2018). 32. Ropo, M., Kokko, K. & Vitos, L. Assessing the Perdew-Burke-Ernzerhof exchange-correlation density functional revised for metallic bulk and surface systems. Phys. Rev. B 77, 195445 (2008). 33. Shibata, M., Yoshida, K. & Furuya, N. Electrochemical synthesis of urea on reduction of carbon dioxide with nitrate and nitrite ions using Cu-loaded gasdiffusion electrode. J. Electroanal. Chem. 387, 143–145 (1995). 34. Wan, H. et al. Electrochemical synthesis of urea: co-reduction of nitric oxide and carbon monoxide. ACS Catal. 13, 1926–1933 (2023). 35. Calle-Vallejo, F., Martínez, J. I., García-Lastra, J. M., Mogensen, M. & Rossmeisl, J. Trends in stability of perovskite oxides. Angew. Chem. Int. Ed. 49, 7699–7701 (2010). 36. Sargeant, E., Illas, F., Rodríguez, P. & Calle-Vallejo, F. Importance of the gasphase error correction for O2 when using DFT to model the oxygen reduction and evolution reactions. J. Electroanal. Chem. 896, 115178 (2021). 37. Calle-Vallejo, F., Huang, M., Henry, J. B., Koper, M. T. M. & Bandarenka, A. S. Theoretical design and experimental implementation of Ag/Au electrodes for the electrochemical reduction of nitrate. Phys. Chem. Chem. Phys. 15, 3196–3202 (2013). 38. Urrego-Ortiz, R., Builes, S. & Calle-Vallejo, F. Fast correction of errors in the DFT-calculated energies of gaseous nitrogen-containing species. ChemCatChem 13, 2508–2516 (2021). 39. Urrego-Ortiz, R., Builes, S. & Calle-Vallejo, F. Impact of intrinsic density functional theory errors on the predictive power of nitrogen cycle electrocatalysis models. ACS Catal. 12, 4784–4791 (2022). 40. Granda-Marulanda, L. P. et al. A semiempirical method to detect and correct DFT-based gas-phase errors and its application in electrocatalysis. ACS Catal. 10, 6900–6907 (2020). 41. Peterson, A. A., Abild-Pedersen, F., Studt, F., Rossmeisl, J. & Nørskov, J. K. How copper catalyzes the electroreduction of carbon dioxide into hydrocarbon fuels. Energy Environ. Sci. 3, 1311–1315 (2010). 42. Studt, F., Abild-Pedersen, F., Varley, J. B. & Nørskov, J. K. CO and CO2 hydrogenation to methanol calculated using the BEEF-VdW functional. Catal. Lett. 143,71–73 (2013). 43. Christensen, R., Hansen, H. A. & Vegge, T. Identifying systematic DFT errors in catalytic reactions. Catal. Sci. Technol. 5, 4946–4949 (2015). 44. Sargeant, E., Illas, F., Rodríguez, P. & Calle-Vallejo, F. On the shifting peak of volcano plots for oxygen reduction and evolution. Electrochim. Acta 426, 140799 (2022). 45. Almeida, M. O., Kolb, M. J., Lanza, M. R. V., Illas, F. & Calle‐Vallejo, F. Gas‐ phase errors affect DFT‐based electrocatalysis models of oxygen reduction to hydrogen peroxide. ChemElectroChem 9, e20220021 (2022). 46. Della Sala, F., Fabiano, E. & Constantin, L. A. Kinetic-energy-density dependent semilocal exchange-correlation functionals. Int. J. Quant. Chem. 116, 1641–1694 (2016). 47. Schmidt, T. & Kümmel, S. Oneand many-electron self-interaction error in local and global hybrid functionals. Phys. Rev. B 93, 165120 (2016). 48. Tao, J., Perdew, J. P., Staroverov, V. N. & Scuseria, G. E. Climbing the density functional ladder: nonempirical meta–generalized gradient approximation designed for molecules and solids. Phys. Rev. Lett. 91, 146401 (2003). 49. Curtiss, L. A., Raghavachari, K., Redfern, P. C. & Pople, J. A. Assessment of gaussian-2 and density functional theories for the computation of enthalpies of formation. J. Chem. Phys. 106, 1063–1079 (1997). 50. Curtiss, L. A., Raghavachari, K., Redfern, P. C. & Pople, J. A. Investigation of the Use of B3LYP zero-point energies and geometries in the calculation of enthalpies of formation. Chem. Phys. Lett. 270, 419–426 (1997). 51. Paier, J., Hirschl, R., Marsman, M. & Kresse, G. The Perdew–Burke–Ernzerh of exchange-correlation functional applied to the G2-1 test set using a planewave basis set. J. Chem. Phys. 122, 234102 (2005). 52. Heyd, J., Scuseria, G. E. & Ernzerhof, M. Hybrid functionals based on a screened coulomb potential. J. Chem. Phys. 118, 8207–8215 (2003). 53. Urrego-Ortiz, R., Builes, S. & Calle-Vallejo, F. Automated versus chemically intuitive deconvolution of density functional theory (DFT)-based gas-phase errors in nitrogen compounds. Ind. Eng. Chem. Res. 61, 13375–13382 (2022). 54. Kresse, G. & Furthmüller, J. Efficient iterative schemes for ab initio totalenergy calculations using a plane-wave basis set. Phys. Rev. B 54, 11169–11186 (1996). 55. Wang, W. et al. Measurement of the cleavage energy of graphite. Nat. Commun. 6, 7853 (2015). 56. Girifalco, L. A. & Lad, R. A. Energy of cohesion, compressibility, and the potential energy functions of the graphite system. J. Chem. Phys. 25, 693–697 (1956). 57. Zacharia, R., Ulbricht, H. & Hertel, T. Interlayer cohesive energy of graphite from thermal desorption of polyaromatic hydrocarbons. Phys. Rev. B 69, 155406 (2004). 58. Benedict, L. X. et al. Microscopic determination of the interlayer binding energy in graphite. Chem. Phys. Lett. 286, 490–496 (1998). 59. Xia, M., Liang, C., Cheng, Z., Hu, R. & Liu, S. The adhesion energy measured by a stress accumulation-peeling mechanism in the exfoliation of graphite. Phys. Chem. Chem. Phys. 21, 1217–1223 (2019). 60. Perdew, J. P. Jacob’s Ladder of Density Functional Approximations for the Exchange-Correlation Energy. In AIP Conference Proceedings. Vol. 577, pp 1–20 (AIP, 2001). 61. Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 77, 3865–3868 (1996). COMMUNICATIONS CHEMISTRY | https://doi.org/10.1038/s42004-023-00990-7 ARTICLE COMMUNICATIONS CHEMISTRY | (2023) 6:196 | https://doi.org/10.1038/s42004-023-00990-7 | www.nature.com/commschem 9