Impact of Intrinsic Density Functional Theory Errors on the Predictive Power of Nitrogen Cycle Electrocatalysis Models.

Oxidized nitrogen species can pollute both the atmosphere[1−3] and water bodies.[4−6] Their concentrations are worryingly increasing because of anthropogenic activities such as the combustion of fossil fuels and intensive agriculture.[7−11] An alternative to remediate their negative impact is to reduce them into unharmful molecular nitrogen (N2) or valuable ammonia (NH3),[12−14] thereby dynamizing the nitrogen cycle. In principle, electrocatalysis could be used as a green technology for these processes if the necessary energy input comes from renewable sources.[12,15] However, the design of active, selective, and stable catalysts for the reduction of nitrogen oxides is not trivial. In that regard, density functional theory (DFT) calculations could serve as a supplement, support, or guide to experiments.[16−23]

DFT is widely used in computational chemistry for the modeling of solids. Specifically, exchange-correlation functionals at the generalized gradient approximation (GGA) have shown high accuracy with low computational requirements when predicting the ground-state properties of bulk and surface metals.[24,25] However, when predicting gas-phase energetics, the limitations of GGA functionals are well-known (e.g., overbinding energy of N2 and O2)[25−27] and predictions in line with experiments are only expected on the basis of error cancellation, i.e., when similar compounds appear in opposite sides of chemical reactions.[28−30] The inaccuracies may be reduced by the use of meta-GGA functionals, which represent a step up in the hierarchy of exchange-correlation approximations.[31] Because functionals at the meta-GGA level take into account the kinetic energy density of the Kohn–Sham orbitals, they are supposedly better than GGAs for molecules, while metals are still accurately described.[32,33]

Gas-phase errors are problematic in heterogeneous catalysis, where an accurate description of the gas phase is paramount for adsorption and desorption steps. Such steps happen each at least once in every catalytic reaction. In spite of their gas-phase errors, GGA functionals are extensively used in catalysis given their low computational requirements. Previous efforts have been devoted to (i) benchmarking their performance for predicting the enthalpies and entropies of adsorption of various systems[34−37] and (ii) combining different functionals to boost their accuracy.[38,39] Considering recent error analysis on nitrogen-containing organic compounds,[40] if DFT at the GGA level is used to model reactions involving nitrogen oxides, it is expected that the calculated energies will entail large errors, in particular for highly oxidized species, such as nitrate and nitrite. Thus, accurately assessing the energetics of reactions such as nitrate reduction or electrochemical nitrogen oxidation remains challenging.

Herein, we show that large errors are encountered in the GGA and meta-GGA formation enthalpies of 11 oxidized nitrogen species in the gas phase. Importantly, the errors scale with the number of oxygens in the structure and the scaling factor is approximately constant for all the functionals studied. This exposes an intrinsic GGA and meta-GGA limitation that must be overcome if accurate predictions are sought after for the modeling of catalytic redox processes among nitrogen-containing species. Furthermore, we show the effects of intrinsic gas-phase errors on adsorption-energy scaling relations and volcano plots for two electrocatalytic reactions and propose an inexpensive scheme to systematically correct such errors.

Computational Details

All the energies were calculated with DFT using the Vienna ab initio simulation package (VASP).[41] The gas-phase calculations were performed for seven exchange-correlation functionals: four GGAs (PBE,[42] PW91,[43] RPBE,[44] and BEEF-vdW[45]), one meta-GGA (TPSS[33]), and two hybrids (PBE0[46] and B3LYP[47]). The adsorption energies on porphyrins were calculated with PBE and RPBE. Molecular representations of the nitrogen-containing compounds studied here are shown in Figure . Besides, Figure S2 provides the skeletal formulas of the oxidized nitrogen species, in which their single and multiple bonds are apparent.

Figure 1

Schematics of the nitrogen species in this work. Purple, red, and pink spheres represent nitrogen, oxygen, and hydrogen atoms, respectively. The skeletal formulas of these molecules, where single and multiple bonds are depicted, can be found in Figure S2.

For metalloporphyrins (see the schematic in Figure S3), spin-unrestricted calculations with and without adsorbates were performed, and the most stable spin state was selected in each case to assess the adsorption energies (Table S4). The computational hydrogen electrode was used to describe the energetics of proton–electron transfers.[48] Further computational details, including the assessment of the free energies of adsorption and a comparison between experimental and computational zero point energies (ZPEs), are provided in section S1 of the Supporting Information. As the experimental and calculated ZPEs are nearly identical, we conclude that the discrepancies in the formation energies stem mainly from the total energies calculated with DFT.

Detection and Correction of the Gas-Phase Errors

The formation energies of nitrogen compounds from their elements in their respective standard states can be calculated from eq .where HNO is an oxidized nitrogen compound. We note that when the compounds do not contain hydrogen (i.e., N2O and NO; see Figure ), x = 0 in eq .

The total errors in the description of the oxidized nitrogen compounds (εH) are determined as the difference between the DFT-calculated and experimental enthalpies of formation (ΔHHDFT and ΔHHexp), as in eq .

Experimental values were taken from thermodynamic tables.[49,50] The total error of HNO encompasses the errors of the reactants and products of eq . Thus, the total error can be estimated from these individual errors as[40]where εH, εN, and εO are the errors of the reactants in eq (H2, N2, and O2) and εH is the gas-phase error of the oxidized nitrogen compound itself, namely, the product of eq . Since H2 is generally well described by DFT, εH ≈ 0. Conversely, the triplet state of O2 is poorly described by GGA functionals,[25] such that εO is typically large.[51] O2 can be swiftly corrected using a semiempirical approach based on the formation energy of H2O.[28,40,48] In addition, εN is usually substantial and can be calculated from the ammonia synthesis reaction, as explained elsewhere[40] and in section S4.

If only the errors in O2 are corrected, the convoluted error of a specific nitrogen compound and that of N2 () can be calculated by combining eqs and 3:

Furthermore, if the errors in O2 and N2 are simultaneously corrected, εH is found to be

Equations and 5 can be used to progressively isolate the errors of all oxidized nitrogen compounds. Figure shows these errors as functions of the number of oxygen atoms in the molecule (nO). Indeed, Figure a (O2 is corrected) and b (O2 and N2 are corrected) shows that such errors are large for GGA and meta-GGA functionals, with values as large as −3.0 eV. More importantly, the errors are linear functions of the number of oxygen atoms in the molecules and the linear trends have, on average, slopes of −0.5 eV/O atom for the GGA and meta-GGA functionals (see the specific values in Table S2). This implies the following:

Figure 2

DFT errors in the formation enthalpies of nitrogen compounds as a function of their number of oxygen atoms (nO). Circles (●) are the calculated data points for the GGAs, squares (■) for TPSS, and triangles (▲) for hybrids. Least-squares linear fits are shown as continuous lines. (a) Errors obtained with the DFT calculations and the corrected energies of O2 using eq (Table S2). (b) Errors obtained after correcting O2 and N2, which correspond to the isolated errors of each nitrogen-containing molecule in eq (Table S3). (c) Residual errors (eq ) after correcting GGA and meta-GGA functionals using eq (Table S11). In all panels the hybrids are used as a benchmark and were not corrected at all (Table S2). The MAEs and MAX for panels a and b are given in section S2. For each value of nO, the species are as follows: 1, N2O, NO, HNO; 2, cis-N2O2, HNO2, NO2; 3, N2O3, HNO3, NO3; 4, N2O4; 5, N2O5.

Progressively adding oxygen atoms to a nitrogen-containing molecule increases the magnitude of the DFT errors by roughly 0.5 eV each time, which is too large for accurate predictions of reaction energies and associated properties such as equilibrium potentials.

The errors are intrinsic; that is, they are due to the generalized gradient approximation. As such, they cannot be avoided by switching to other GGA or meta-GGA functionals. Conversely, the hybrid functionals PBE0 and B3LYP were not corrected at all because their trends in Figure display nearly flat slopes of −0.06 and 0.04 eV/O atom, respectively. Besides, they have mean absolute errors (MAEs) of 0.16 and 0.10 eV. This is consistent with hybrid functionals being generally able to reproduce the experimental energetics of small molecules more closely than GGAs.[27,46,52,53]

DFT-based modeling of redox processes among nitrogen-containing compounds in Figure entails sizable errors, in particular when there are large differences in the oxidation states of the reactants and products (e.g., nitrate reduction to N2).

Because the errors are systematic, a model can be made that simultaneously corrects all errors based on nO. One such method is detailed in the next paragraphs. We note that the dependence of the DFT errors on nO can be rationalized by the presence of multiple (i.e., double or triple) bonds in HNO, since it is known that DFT-GGAs often fail to accurately describe molecules with such bonds.[25,36]Figure S2 shows the skeletal structures of the molecules in Figure , in which single and multiple bonds are apparent. As all of them have unsaturated bonds and, in several cases, single and multiple bonds are intercalated, resonant structures are possible, which likely induce the large errors observed.[54] This is in line with previous works showing that compounds comprising multiple bonds, such as nitrates or carboxylic acids, display large errors.[30,40]

When the linear trends in Figure a are used to correct the intrinsic errors in the formation enthalpies of the nitrogen compounds of the GGA and meta-GGA functionals, the averages of the MAEs and maximum absolute errors (MAX) are 0.18 and 0.40 eV (Table S12). Similarly, if the linear trends in Figure b are used to correct the calculated formation enthalpies, large errors are also obtained, with the averages of the MAEs and MAX being 0.18 and 0.38 eV (Table S13). The considerably lower errors in Figure c (average MAE of 0.04 eV and average MAX of 0.12 eV) are obtained by splitting the nitrogen-containing molecules into the three groups shown in Figure : (i) dinitrogen-containing species (N2O), (ii) mononitrogen-containing species (NO), and (iii) hydrogenated species (HNO). We note that similar categories have previously been used to rationalize energetic and structural differences of oxidized nitrogen species.[55] For each of these three groups, nO is still linearly related to the errors and can be used to correct the DFT-calculated enthalpy (Figure S1). In this order of ideas, the corrected formation enthalpies are given by eq :where ΔHcorr is the corrected enthalpy of the oxidized nitrogen species and m and b are, respectively, the slope and intercept of the regression line of group i = N2O, NO, and HNO. The functional-dependent values of m and b are reported in Table . Regardless of the functional, the NO group has steeper slopes compared to the N2O and HNO groups. In fact, the average slopes for the three groups are −0.67 (NO), −0.43 (N2O), and −0.45 eV/O atom (HNO). Further details of the fitting procedure appear in section S2. In analogy to eq , we calculate the residual errors (ε′) as

Table 1

Parameters to Correct the Formation Enthalpies of the Nitrogen-Containing Species for Each GGA and meta-GGA Functional Studied Using Eq a

Parameter	PBE	PW91	RPBE	BEEF-vdW	TPSS
mN2Ox	–0.42	–0.40	–0.41	–0.46	–0.49
bN2Ox	0.00	0.10	–0.46	–0.67	–0.61
mHNOx	–0.45	–0.43	–0.42	–0.46	–0.50
bHNOx	0.39	0.41	0.10	0.01	0.04
mNOx	–0.67	–0.66	–0.65	–0.67	–0.68
bNOx	0.59	0.69	0.23	0.09	0.05

The slopes (m) are in eV/O atom, and the intercepts (b) are in eV.

We take N2O5 calculated with PBE to illustrate the use of the corrections from Table and eq . In this case, the experimental value is ΔHNexp = 0.12 eV and the DFT-calculated formation enthalpy after the O2 correction is . This nitrogen oxide belongs to the group of dinitrogen oxides (N2O); thus, y/2 = 1 and εNPBE = 0.34 eV.[40] Hence, we have the following:. N2O5 belongs to the group of dinitrogen-containing species (N2O) and contains five oxygen atoms; thus, nO = 5. From Table for N2O, m = −0.42 eV/O and b = 0.00 eV. Hence, ΔHNcorr, PBE = −1.93–5·(−0.42) – (0.00) = 0.16 eV, which deviates from experiments by 0.04 eV (εN′ = 0.04 eV). Thus, after applying the corrections, the error changes from εN = −2.38 eV to εN′ = 0.04 eV.

Figure shows for all functionals under study the MAEs and MAX upon subsequently applying these corrections to the nitrogen species in Figure . The final MAEs and MAX of the corrected GGA and meta-GGA functionals are smaller than those of the hybrids and in all cases close to chemical accuracy (1 kcal/mol, red lines in Figure ). In Figure , correcting the error in N2 does not necessarily improve the gas-phase errors. Indeed, for PBE and PW91, the errors are lowered after correcting N2, but the values for RPBE, BEEF-vdW, and TPSS increase. This behavior is not random but depends on εN: if it has the same sign as the errors of the oxidized nitrogen species, εN cancels out a portion of those (eq ).[40] In addition, the small change for RPBE stems from its small εN of −0.05 eV.

Figure 3

(a) Mean and (b) maximum absolute errors (MAE and MAX) for selected functionals after correcting O2 (black), O2 and N2 (gray), and after applying the correction method based on Table and eq (white) to the nitrogen species in Figure . The results for hybrid functionals (B3LYP and PBE0) are used as a benchmark, so that they were corrected in neither panel a nor panel b. The red line represents chemical accuracy (1 kcal/mol).

Impact on Catalysis

The errors in the previous section are calculated only for gaseous compounds and not for clean or adsorbate-covered active sites. In the following, we will assume that the errors of the active sites with and without adsorbates are similar. This was shown to be a good approximation for the modeling of CO2 electroreduction to CO on Cu, Ag, and Au electrodes after applying gas-phase corrections.[30] However, we cannot discard the idea that significant errors might in some cases subsist after correcting the gas phase.[56]

The importance of the proposed gas-phase corrections for the modeling of catalytic processes within the N cycle is apparent when the adsorption-energy scaling relations[57−59] among nitrogen oxides are considered. Figure provides the free energies of adsorption of NO3 and NO2 as a function of that of NO on the metal center of six porphyrins with MN4 sites (M: Ti, V, Cr, Mn, Fe, and Co) using RPBE.

Figure 4

Adsorption energies of (a) NO3 and (b) NO2 as a function of those of NO on M-porphyrins (M = Ti, V, Cr, Mn, Fe, and Co) calculated with RPBE. Dashed lines and open circles correspond to uncorrected DFT calculations, whereas solid lines and full circles correspond to the results upon gas-phase corrections. The brown arrows show the magnitude and direction of the corrections εH for V porphyrin, which are identical for the rest of the materials in the trends.

Since all the species involved in the scaling relations have an associated gas-phase error, once their energies are corrected, each point is vectorially displaced in the plot. Figure illustrates this for the particular case of a V porphyrin. In Figure a the NO correction displaces the data point to the left by 0.41 eV (εNORPBE = −0.41 eV; see Table S3), while the NO3 correction displaces the data point downward by 1.72 eV (εNORPBE = −1.72 eV; Table S3), resulting in a net diagonal displacement of 1.77 eV. In Figure b the NO correction is identical and that of NO2 displaces the data point downward by 1.12 eV (εNORPBE = −1.12 eV; Table S3), resulting in a net diagonal displacement of 1.19 eV. Since each point in the scaling relation is shifted by a constant amount after applying the corrections, the slopes of the scaling relations remain constant, but the intercepts change. For instance, in Figure b the intercept is initially 0.88 eV, and the slope of the uncorrected scaling relation is 1.48 eV. The errors in NO2 and NO are −1.12 and −0.41 eV, so that upon corrections one gets an offset of (0.88 – 1.12 + 1.48·0.41) = 0.37 eV.

Scaling relations are extensively used to build the so-called Sabatier volcano plots.[58−60] Those activity plots help find the adsorption energies of key intermediates that ensure optimal catalysis. Because gas-phase corrections modify the offsets of scaling relations (as in Figure ), the volcano plots based on them are appreciably modified as well. This is exemplified for the electrochemical ammonia synthesis reaction (N2 + 6H+ + 6e– → 2NH3) on metalloporphyrins in Figure calculated with PBE.

Figure 5

Volcano plot for the electrochemical ammonia synthesis for metalloporphyrins using PBE. We provide a wide-range analysis in panel a, while a focus into the dashed region around the volcano tops is shown in panel b. Red lines and open circles correspond to the uncorrected DFT calculations. Green lines and solid circles represent the results upon correcting the gas-phase errors of N2. The red/green dashed lines are the equilibrium potential before/after correcting the N2 errors. The arrows in panel b indicate that to reduce the limiting potential, *NH2 binding on V porphyrin has to be weakened or strengthened depending on the inclusion or exclusion of gas-phase corrections.

First of all, in Figure a the equilibrium potential of electrochemical ammonia synthesis is presented before (red dashed line) and after (green dashed line) correcting the N2 error. Before correcting the equilibrium potential, PBE yields 0.113 V vs RHE. Once the proposed corrections are applied, the equilibrium potential is 0.057 V vs RHE, in agreement with experiments. Admittedly, the error in the equilibrium potentials is not large, but it is amplified by a factor of 6 when assessing the corresponding reaction energies. For PBE this results in a deviation with respect to experiments of the reaction energy as large as 0.34 eV.

As each side of the volcano corresponds to a different electrochemical step, the gas-phase corrections are different. Specifically, the potential on the right leg of the volcano is typically limited by[23] N2+* + H+ + e– → *N2H, which involves the error in N2. In turn, the potential on the left leg is usually limited by *NH2 + H++e → NH3+*, which means that the values on that side need not be corrected, as NH3 is generally well described by DFT.[25] Since the legs of the Sabatier volcanoes are based on scaling relations, their slopes are identical before and after correcting the gas-phase errors, yet their offsets vary. This causes the following:

(i) A change in the limiting potential (U in V vs RHE; see Table S15). For instance, the limiting potentials of the corrected porphyrins on the right leg of the volcano in Figure a (in green) are shifted downward by 0.34 eV with respect to the uncorrected data points (in red), because of the error in N2. However, Ti porphyrins experience no change in U. In turn, V porphyrins experience an intermediate downward shift of 0.15 eV stemming from the switch of potential-limiting steps (Figure b).

(ii) A shift in the location of the volcano apex. Such a shift can change activity orderings, the specific volcano leg a material belongs to, and the energy difference between the calculated data points and the apex. This is apparent in Figure b for V porphyrin, which switches from the left to the right of the volcano when the gas phase is corrected (see section S5). Before the corrections, V porphyrin is below the volcano apex by 0.09 V. Upon the corrections, it is below by 0.03 V. In terms of ΔGNH, before the corrections V porphyrin is 0.09 eV to the left of the apex. Afterward, it is 0.12 eV to the right. This means that the guidelines for optimizing this porphyrin derived from the two volcanoes are the exact opposite. For comparison, an analogous analysis is provided in Figure S4 for RPBE. We conclude that large shifts in the top of the volcano are normally associated with large gas-phase errors. We emphasize that knowing the precise location of a material with respect to the volcano apex is crucial for its optimization. This has been profusely illustrated for Pt-based catalysts for O2 reduction, where Pt(111) is located 0.1 eV to the left of the top in terms of *OH binding energy,[61,62] such that catalysts are engineered to bind *OH more weakly than Pt(111) by no more than 0.1 eV.

To further illustrate the effect of gas-phase errors in electrocatalysis of the nitrogen cycle, we used RPBE with and without gas-phase corrections to model nitric oxide reduction to hydroxylamine (NO(g) + 4H+ + 3e– → NH3OH+(aq)) on porphyrin catalysts. Hydroxylamine is a value-added chemical with numerous uses in industry, such that its electrochemical production from nitrate or NO is an economically appealing way of balancing the N cycle.[63,64] Generally, we find that the lowest-energy pathway on metalloporphyrins is NO → *NHO → *ONH2 → NH3OH+ (see Table S16).

The energetics of protonated hydroxylamine (NH3OH+(aq)) was calculated on the basis of its acid–base equilibrium with neutral hydroxylamine (NH2OH, pKa = 7.68),[65] the formation free energy of which in the gas phase is −0.03 eV.[66] Besides, we found a gas-phase error in NH2OH of −0.13 eV. As shown in Figure , the volcano plot for this reaction has the usual right (weak binding) and left (strong binding) regions and an intermediate binding region. In Figure , we observe the following:

Figure 6

Volcano plot for NO reduction to hydroxylamine on metalloporphyrins calculated with (a) RPBE and (b) RPBE and gas-phase corrections. The red and green dashed lines represent the calculated equilibrium potential, while the black dashed line represents the experimental value. Materials on the left are generally limited by *ONH2 hydrogenation, those on the right by NO hydrogenation to *NHO, and those in the middle by *NHO hydrogenation to *ONH2; see Table S16.

(i) The equilibrium potentials and reaction energies change upon applying gas-phase corrections: 0.36 and 0.45 V vs RHE for uncorrected and corrected RPBE, while the experimental value is 0.46 V vs RHE. The respective reaction energies are −1.08, −1.36, and −1.38 eV.

(ii) The ordering of catalytic activities changes upon applying gas-phase corrections. For uncorrected RPBE it is Cr > Co > Mn > Fe > V > Ni > Cu > Ti. Conversely, the ordering for corrected RPBE is Co > Cr > Mn > Fe > Ni > Cu > V > Ti.

(iii) The potential-limiting steps of some materials change upon applying gas-phase corrections. This is the case of the most active porphyrins, namely those of Cr (step 2 vs 3 in uncorrected and corrected RPBE) and Co (step 1 vs 2). Such changes are worth noting, as the recipes for optimizing those materials change significantly (as shown before in Figure b).

(iv) Quantitative differences are found for the limiting potentials of the electrocatalysts under study, as shown in Tables S17–S19. The differences are linked to the gas-phase corrections of NO and NH2OH (−0.41 and −0.13 eV), which is apparent when the potential-limiting steps are the same with and without gas-phase corrections.

(v) The length of the intermediate binding region changes upon gas-phase corrections. In terms of ΔGONH, such region has a length of 0.25 and 0.43 eV without and with corrections, respectively.

Concluding Remarks

Using DFT to predict the gas-phase formation enthalpies of oxidized nitrogen species often results in large negative errors when using GGA and meta-GGA functionals. The DFT errors scale approximately linearly with the number of oxygen atoms in the compounds, with slopes for GGA and meta-GGA functionals close to −0.5 eV/O atom. This is considerably steeper than those found for hybrid functionals (−0.06 and 0.04 eV/O atom for PBE0 and B3LYP). The similar slopes among GGA and meta-GGA functionals indicate that DFT predictions worsen progressively when adding O atoms to the structure and suggest that such errors are probably unavoidable at those rungs of Jacob’s ladder of density functional approximations.[31]

If the data set is subdivided into dinitrogen-containing species (N2O), mononitrogen-containing species (NO), and hydrogenated species (HNO), the resulting formation enthalpies can be swiftly corrected with high accuracy. In fact, the MAEs of the GGA and meta-GGA functionals are reduced, on average, from 1.23 to 0.04 eV.

Furthermore, gas-phase errors significantly alter adsorption-energy scaling relations and the volcano plots built upon them. We exemplified that for electrochemical ammonia synthesis and NO reduction to hydroxylamine. The magnitude and direction of the displacements depend on the separate gas-phase errors and the slope of the scaling relations. We noticed changes in (i) the equilibrium potentials and reaction energies, (ii) the location of the volcano peaks, (iii) the location of catalysts on the regions of the volcano, (iv) the predicted limiting potentials, and (v) the catalytic activity orderings. The changes have a direct connection with the magnitudes of the gas-phase errors.

We hope that, as the electrocatalysis of the N cycle regains more and more attention, computational chemists will become increasingly aware of the fact that GGA and meta-GGA functionals have intrinsic gas-phase errors that may impair their catalytic activity, selectivity, and stability predictions. Finally, we note that further experimental and computational efforts are necessary to detect and correct possible errors in the active sites with and without adsorbates.