Automated versus Chemically Intuitive Deconvolution of Density Functional Theory (DFT)-Based Gas-Phase Errors in Nitrogen Compounds.

Catalysis models involving metal surfaces and gases are regularly based on density functional theory (DFT) calculations at the generalized gradient approximation (GGA). Such models may have large errors in view of the poor DFT-GGA description of gas-phase molecules with multiple bonds. Here, we analyze three correction schemes for the PBE-calculated Gibbs energies of formation of 13 nitrogen compounds. The first scheme is sequential and based on chemical intuition, the second one is an automated optimization based on chemical bonds, and the third one is an automated optimization that capitalizes on the errors found by the first scheme. The mean and maximum absolute errors are brought down close to chemical accuracy by the third approach by correcting the inaccuracies in the NNO and ONO backbones and those in N-O and N-N bonds. This work shows that chemical intuition and automated optimization can be combined to swiftly enhance the predictiveness of DFT-GGA calculations of gases.

Introduction

Because of its wide range of oxidation states from −3 to +5, nitrogen forms a wide and diverse group of compounds when combined with hydrogen and oxygen, including oxides, hydrides, radicals, ions, and acids.[1] All those compounds are part of the nitrogen cycle and are relevant in aquatic and terrestrial systems, atmospheric chemistry, and chemical industries.[2−6] The reactions connecting these compounds have gained interest in the scientific community, because of their industrial uses, the adverse effect some of them have on human health, their role in climate change, and the colossal imbalance of the nitrogen cycle as a result of human activities.[1,7−11]

The experimental assessment of the chemical properties of nitrogen-containing species is far from straightforward, given their instability/reactivity and the complex reaction networks they form.[12−14] Hence, density functional theory (DFT) calculations have been widely used to predict structural and thermochemical properties of this family of compounds.[15−18] To achieve fair predictions of molecular systems, DFT calculations frequently rely on hybrid exchange-correlation functionals, such as B3LYP,[19] which are computationally expensive and not advisible for systems with delocalized electrons.[20,21] This is certainly restrictive for studies in areas such as heterogeneous electrocatalysis, where the systems often involve gases and liquids in contact with conductive solids.[22]

Exchange-correlation functionals based on the generalized gradient approximation (GGA) are widely used to study catalytic systems in view of their affordable computational requirements and reasonable predictions of metal bulk and surface properties.[23−26] Although the limitations of GGAs in describing the energetics of molecules are well-known,[27−29] semiempirical and fully computational schemes can be devised to rapidly correct them.[30−34] Correction schemes generally seek to identify specific groups of atoms or functional groups that systematically contribute to the total errors of different molecules. While the errors can mostly be attributed to poor descriptions of the exchange contribution to the total DFT energy in molecules with multiple bonds and/or strongly interacting lone pairs,[27,35] the detection and quantification of systematic errors enables otherwise modest GGA functionals to produce accurate yet inexpensive predictions.

In this paper, we analyze three approaches to determine the error contributions in the formation energies of 13 compounds containing N, O and H, using the GGA-PBE exchange-correlation functional.[36] The first approach is sequential and based on chemical intuition, the second approach is automated and based on chemical bonds, and the third approach is also automated and capitalizes on the findings of the first one. The sequential approach finds the errors in a stepwise fashion, while the automated approaches minimize error functions with free parameters. Apart from swiftly bringing the PBE-calculated Gibbs energies close to the experimental values, our results show how chemical intuition can be used as an initial step to identify errors that automated schemes can further minimize.

Computational Methods

Ball-and-stick representations of the nitrogen-containing compounds studied here (NH2OH, NO, HNO, NO2, NO3, trans-HNO2, cis-HNO2, HNO3, N2O, cis-N2O2, N2O3, N2O4, N2O5) are shown in Figure S1 in the Supporting Information. The DFT calculations were performed using the Vienna ab initio simulation package (VASP),[37] the Perdew–Burke–Ernzerhof (PBE) exchange-correlation functional,[36] and the projector augmented-wave (PAW) method.[38] The plane-wave cutoff for all the calculations was 450 eV, shown previously[30,39,40] and verified in Figure S2 in the Supporting Information for N2O formation to provide converged reaction energies. Gaussian smearing with kBT = 10–3 eV was used and all energies were extrapolated to 0 K. During the structural optimization of the molecules, carried out using the conjugate gradient algorithm, all atoms were allowed to relax in all directions until the maximal atomic forces were equal to or smaller than 0.01 eV/Å. The molecules were simulated in large boxes in which the distance between periodic images was at least 12 Å. Accordingly, we only considered the Γ-point for the k-point sampling of the calculations. Spin-unrestricted calculations were performed for O2, NO, NO2, NO3, and cis-N2O2. Note that NO, NO2, and NO3 are neutral free radicals, not anions or cations. The error optimizations were formulated in GAMS using the ANTIGONE[41] global optimization solver on the NEOS Server.[42,43]

The formation of a generic nitrogen compound HNO from its elements in their respective standard states is defined as

If the molecule does not contain hydrogen (e.g., NO3), x = 0 in eq . Likewise, if the molecule does not contain oxygen (e.g., NH3), z = 0 in eq . Following previous works,[30,32,39,40] the total error in the DFT description of a generic nitrogen compound (εHT) is defined as the difference between the DFT-calculated and the experimental energies of formation, see eq . In this case, we will make the analysis in terms of Gibbs energies of formation (ΔfGHDFT and ΔfGHexp). However, we note that the analysis can be made in terms of enthalpies of formation and the results would be identical, because the total entropies of the molecules are usually taken from tabulated experimental data.[44,45]

The Gibbs energies of formation were approximated by means of DFT as follows:where ΔfEDFT is the formation energy calculated with DFT total energies, ΔfZPE the zero-point energy change calculated using DFT within the harmonic oscillator approximation, and TΔfS the entropy change at T = 298.15 K, taken from thermodynamic tables.[44,45] We did not incorporate heat capacity contributions to the formation energies in eq , because their energy change has been shown to be small in the range of 0 to 298.15 K[31,46] (see further details in section S6 in the Supporting Information). We note that previous works showed that the differences between experimental and calculated ΔfZPE are negligible for various HNO compounds,[34] such that the errors can be entirely assigned to ΔfEDFT. The experimental Gibbs energies used in eq to compute the errors are also taken from thermodynamic tables.[44,45]

As shown in eq , the total error in eq for HNO is the difference in the errors of the products and reactants as given by eq :

In eq , εH, εN, and εO are the respective individual errors of H2, N2, and O2, and εHNO is the gas-phase error of the generic nitrogen compound. A usual simplification is εH ≈ 0, because H2 is generally well described by DFT.[27] In contrast, εN and εO are generally large[27,47] and are respectively assessed based on the reactions in which N2 is combined with H2 to produce NH3,[32] and O2 is combined with H2 to produce H2O,[39,48]In those reactions, namely, ammonia synthesis and water formation, εN and εO can be isolated because both ammonia and water are generally well described by DFT, as they only contain N–H and O–H single bonds.[27] Once the errors in O2 and N2 are corrected, εHNO can be calculated by combining eqs and 4:

As will be shown in the subsequent section, there are ways of estimating εHNO based on the bonds and/or groups of atoms present in HNO. The accuracy of the εHNO estimates will dictate the magnitude of the final errors (calculated with an updated version of eq ), see below. In the following, we assume that εHNO is the sum of the errors due to the bonds and/or groups of atoms in HNO. While this is usually a fair assumption, previous works[32] showed that if a large functional group is present more than once in a small compound, additional intramolecular interactions might appear that change the magnitude of gas-phase corrections. In any case, it is possible to use the estimates of εHNO to correct the DFT-calculated Gibbs energy of formation as

Results and Discussion

Once the errors have been isolated for each compound using eq , we employ a sequential correction approach (herein referred to as sequential) based on the complexity of the molecules (see the flowchart of the method in Figure ). By complexity, we mean an increasing number of bonds and/or groups of atoms in the molecules. Initially, the simplest possible molecules are analyzed, and then increasingly large molecules are considered, to detect the possible bonds and/or functional groups responsible for the errors in their structures. In that order of ideas, the first compound we analyzed was hydroxylamine (NH2OH), as it only features single N–H, N–O and O–H bonds. The difference between the calculated and experimental formation energies is −0.15 eV, from which we conclude that single N–O bonds ought to be corrected by that much.

Figure 1

Flowchart of the sequential method used to correct the gas-phase formation energies of HNO. Increasingly complex molecules are analyzed in every step until the entire dataset is corrected.

Next, we analyzed the N=O double bond present in nitric oxide (NO) and nitroxyl (HNO), the formation energies of which differ from experiments by 0.07 and 0.04 eV. As these errors are smaller than 0.10 eV, we opt not to correct them, although an additional average correction is an option if results with higher accuracy were needed. The next compound in the list is nitrogen dioxide (NO2), which displays a large error with respect to experiments of −0.79 eV. From this we conclude that the presence of an ONO backbone in a molecule induces an error of −0.79 eV in its calculated formation energy. Analogously, we note that the OCO backbone has been previously identified to induce appreciable errors in organic and inorganic C-containing compounds.[30,31] The ONO backbone error together with the N–O single bond error allow us to correct the formation energies of nitrogen trioxide (NO3), trans and cis nitrous acid (trans-HNO2, cis-HNO2), and nitric acid (HNO3), which initially differ from experiments by as much as 1.38, 0.52, 0.52, and 0.94 eV. Upon the corrections, the residual errors are 0.19, 0.02, 0.02, and 0.00 eV.

The next compound in the list is nitrous oxide (N2O), the DFT formation energy of which initially differs from experiments by −0.49 eV. Therefore, we conclude that molecules with an NNO backbone have an error associated with it as large as −0.49 eV. In view of the lack of information, we assumed that errors in N–N bonds are one-half of those in the NNO backbone and the validity of this assumption will be asserted later in this work. Correcting the NNO, ONO, N–O and N–N errors, we are able to considerably lower the total errors in the formation energies from 0.83 to 0.10 eV for cis-N2O2, from 1.21 to 0.07 eV for N2O3, from 1.80 to 0.02 eV for N2O4, and from 1.98 to 0.11 eV for N2O5.

Overall, the mean and maximum absolute errors (MAE and MAX) are initially 0.83 and 1.98 eV. Once the N–O, N–N, ONO, and NNO errors have been corrected, the resulting MAE and MAX are drastically reduced to 0.07 and 0.19 eV. The initial and final errors for all compounds under study are provided in Figure .

Figure 2

Initial (red) and final (green) errors in the nitrogen compounds under study. The final errors are obtained upon applying a sequential method that identifies N–OH and N–N bonds and ONO and NNO groups as sources of error. The mean and maximum absolute errors (MAE and MAX) are provided before and after the corrections in the inset. The experimental and final Gibbs energies are given in Table .

Table 1

Formation Energies of Nitrogen Compoundsa

	Formation Energies (eV)
species	ΔfGHxNyOzexp	ΔfGHxNyOzONC	ΔfGHxNyOzseq	ΔfGHxNyOzAO1	ΔfGHxNyOzAO2
NH2OH	0.04	–0.11	0.04	0.13	–0.01
NO	0.91	0.83	0.83	1.12	0.83
HNO	1.16	1.12	1.12	1.41	1.12
NO2	0.53	–0.26	0.53	0.44	0.61
NO3	1.20	–0.18	1.01	0.95	1.13
trans-HNO2	–0.46	–0.98	–0.44	–0.46	–0.46
cis-HNO2	–0.43	–0.95	–0.41	–0.43	–0.43
HNO3	–0.76	–1.70	–0.76	–0.76	–0.74
N2O	1.07	0.59	1.07	1.14	1.01
cis-N2O2	2.22	1.39	2.11	2.22	2.16
N2O3	1.48	0.27	1.54	1.52	1.55
N2O4	1.03	–0.77	1.05	0.90	1.03
N2O5	1.21	–0.77	1.11	1.47	1.14

Second column (energies denoted with a superscript “exp”) shows experimental formation energies; the third column (energies denoted with a superscript “ONC”) contains the DFT-calculated formation energies with O2 and N2 corrections. The fourth, fifth, and sixth columns give the final formation energies upon applying the sequential method in Figures and 2 (energies denoted with a superscript “seq”) and automated optimizations 1 and 2 in Figures and 4, respectively (energies denoted with superscripts “AO1” and “AO2”).

A different approach consists of treating the problem as a mathematical optimization in which all errors are simultaneously lowered (hereon referred to as automated optimization 1, AO1). This strategy requires no chemical intuition to hierarchize the chemical compounds in increasing order of complexity, as in the sequential method. AO1 is a multiobjective optimization problem, where the MAE and MAX are minimized simultaneously (for further details see section S4 in the Supporting Information). The adjustable parameters for the optimization can be the errors in single N–O bonds, double N–O bonds, N–N bonds, and O–H bonds. A matrix can be built that decomposes every molecule into these bonds, such that the total error is a sum of all those contributions, see the representation of the 13 nitrogen compounds under study in section S3 in the Supporting Information.

In this case, there is no Pareto front as the MAE and MAX are minimized at the same point. The final MAE and MAX after AO1 are 0.11 and 0.26 eV. The initial and final errors upon this optimization are provided in Figure . Although this procedure substantially lowers the initial MAE and MAX, the sequential method in Figures and 2 performs better (final MAEs: 0.07 vs 0.11 eV; final MAXs: 0.19 vs 0.26 eV).

Figure 3

Initial (red) and final (orange) errors in the nitrogen compounds under study for AO1. The final errors are obtained upon minimizing simultaneously the MAE and MAX using as adjustable parameters the errors in N–O (single and double), N–N and O–H bonds. Inset: MAE and MAX before and after the corrections. The experimental and final Gibbs energies are given in Table .

Similarly, one can capitalize on the errors found by the sequential method by using them as adjustable parameters for another optimization (hereon referred to as automated optimization 2, AO2), see sections S3 and S4. For this multiobjective optimization, the minimum distance selection method[49] was used to find the most feasible point among the Pareto front, see section S4. The knee point or most satisfactory solution inside the feasible space corresponds to a MAE of 0.05 eV and a MAX of 0.08 eV. The initial and final errors obtained after this optimization are provided in Figure . The final MAE and MAX are visibly lower in Figure (0.05 and 0.08 eV), compared to Figures (0.07 and 0.19 eV) and 3 (0.11 and 0.26 eV), and a compromise between the magnitude of the MAX and MAE is attained. Section S5 in the Supporting Information also shows that AO2 is more accurate than a previous method based on the number of oxygen atoms in the molecules.

Figure 4

Initial (red) and final (cyan) errors in the nitrogen compounds under study for AO2. The final errors are obtained from the knee point of the simultaneous minimization of the MAE and MAX using as adjustable parameters the errors in N–O and N–N bonds and those in ONO and NNO groups. Inset: MAE and MAX values before and after the corrections. The experimental and final free energies are given in Table .

In Figures –4 we observe, in broad terms, that the larger the molecules, the larger the errors. In principle, this is because larger molecules have at least one problematic bond and/or group of atoms. However, this trend is not uniform neither in the initial nor the final errors, as the matrix representations of molecules that differ even by one oxygen atom can be rather different. For instance, N2O4 and N2O5 both have two ONO groups in the matrix representation of the sequential method, but the former has a NN bond while the latter has two NO single bonds (see Table S3). Within AO1, N2O4 and N2O5 both have two N=O bonds but the former has a N–N bond and two N–O bonds while the latter has four N–O bonds (see Table S4 in the Supporting Information).

Apart from visualizing the initial and final errors as in Figures –4, it is also convenient to analyze the formation energies to draw further conclusions. All formation energies are provided in Table , and Figure presents a parity plot in which the experimental and calculated (using eq ) formation energies are compared. The fact that the initial (red) data are below the parity line is explained in section S7 in the Supporting Information. The largest initial and final errors appear mostly in the range of 0.80 to 1.65 eV, which contains the lesser stable nitrogen compounds. cis-N2O2, which is the least-stable compound in this study, is an exception, because it displays large initial errors but small final errors. Importantly, the largest final errors in the sequential method and the two automated optimizations correspond to NO3. This suggests that the matrix representations of this compound might be somehow incomplete and/or that it might be advisible to assign a specific error to it if higher accuracy is needed.

Figure 5

Parity plot for the formation energies of the nitrogen compounds under study. Red denotes data without corrections; green represents data corrected by the sequential method (Figure ); orange denotes results from an automated optimization using as free parameters the errors in NO (single and double), NN and OH bonds (AO1, Figure ); cyan represents the results from an automated optimization using as free parameters the errors in N–O and N–N bonds and in the ONO and NNO groups (AO2, Figure ). Inset: region with the largest initial errors (ΔfGexp from 0.80 to 1.65 eV). The gray band covers ±0.20 eV around the parity line.

Table shows the errors found by the three approaches for the bonds and groups of atoms within HNO. The similarities between the values for the sequential method and AO2 are apparent and lead to three observations:

Table 2

Errors Present in Nitrogen Compounds, As Predicted by Three Different Approaches on the Basis of the Bonds and Groups of Atoms Present in the Moleculesa

	Errors (eV)
method	N–H	N–O	N=O	N–N	O–H	NOH	ONO	NNO
AO1	0.00	–0.42	–0.29	–0.26	0.18	–	–	–
AO2	0.00	–	–	–0.07	–	–0.09	–0.87	–0.42
sequential	0.00	–	–	–0.24	–	–0.15	–0.79	–0.49

AO1 and AO2 are shown in Figures and 4, and the sequential method is shown in Figures and 2. Further details appear in sections S3 and S4 in the Supporting Information.

Chemical intuition is able to detect and correct the errors to a great extent.

Averaging among similar compounds leads to even greater accuracy, in particular in this case by lowering the MAX.

Approximating the error in N–N bonds as one-half of the NNO error, as done in the sequential method, is the main difference with respect to AO2.

Based on these observations, we first averaged the errors of the sequential method for NO2 and NO3 to assess the ONO error. Second, N2H4, which should only have an error in its N–N bond, was calculated in a previous work.[32] The error of −0.09 eV reported in that work is close to that found by AO2 (−0.07 eV, see Table ). With these two amends, the MAE and MAX of the sequential method are lowered to 0.06 and 0.13 eV.

Furthermore, some bond errors in AO1 can be combined to approximate the group errors found by the sequential method. For instance, summing the errors of N–O and N=O bonds gives −0.70 eV, which is close to the value of −0.79 eV obtained by the sequential method for the ONO group. Besides, the sum of the N=O and N–N bonds found with AO1 is −0.55 eV, which is also close to the value obtained by the sequential method for the NNO group (−0.49 eV). The same holds for the NOH group: from the addition of the N–O and O–H errors we get −0.23 eV, while the sequential method finds −0.15 eV. Lastly, the N–N errors between the sequential method and AO1 are also rather close (−0.24 and −0.26 eV). However, we stress that it is the small discrepancies between the methods that are ultimately responsible for the different final MAEs and MAXs observed in Figures –4. We also note that section S5 in the Supporting Information verifies that N–H bonds are generally well described and expands on the error in O–H bonds found by AO1.

Before closing the discussion, we note that, when the decomposition of a given molecule into its groups is not unambiguous, running several tests is advisible to determine the best representation. For instance, trans-HNO2 in the sequential method can be thought of having an ONO group and an OH bond. This leads to a residual error of 0.35 eV. In contrast, considering it to be composed of 0.5 ONO groups and a single NO bond, the residual error is 0.00 eV. The latter representation can also be used for cis-HNO2 and extended to HNO3. Moreover, cis-N2O2 can be represented as two NO units linked by a NN bond, leading to a residual error of 0.76 eV. When represented as two NNO units and the double-counting of the NN bond is discounted, the residual error is 0.06 eV.

Impact on Heterogeneous (Electro)catalysis

By correcting gas-phase errors, it is possible to obtain without relying on fortuitous error cancellation accurate reaction energies, equilibrium potentials and, in principle, adsorption energies. This is critical for the models used in computational heterogeneous (electro)catalysis, which are usually based on these properties.

For example, it was shown in a recent work for electrochemical ammonia synthesis and electrochemical nitric oxide reduction to hydroxylamine that gas-phase corrections modify the predicted overpotentials, Sabatier-type volcano plots, and the ordering of catalytic activities among the analyzed materials.[34] In addition, gas-phase corrections have also been shown to improve the prediction of equilibrium and onset potentials for the electroreduction of CO2 to CO,[30] and brought DFT-calculated adsorption energies of CO on several metals closer to experimental values.[50] Finally, the importance of gas-phase corrections has also been illustrated for free-energy diagrams and volcano plots for O2 reduction and evolution[39,51] and H2O2 production.[40]

Furthermore, if a given compound X participates in a catalytic reaction but no experimental data are available for it, one can decompose it in its bonds and/or groups of atoms and anticipate the errors present in its DFT-calculated free energy of formation. However, it is recommendable to make an ensemble of predictions based on different matrix representations and establish a correction range rather than a specific correction.

Finally, we believe that the proposed approaches could be transferred to assess the errors of adsorbates at surfaces. Given their semiempirical nature, this extension would presuppose the availability of accurate experimental adsorption energies.

Conclusions

Herein, we showed that large errors are found when using PBE to assess the free energies of formation of 13 gaseous compounds containing nitrogen, hydrogen, and oxygen. To identify and reduce such errors, we proposed approaches based on chemical intuition and a matrix representation of the molecules. The representations decompose each molecule into the bonds and/or groups of atoms they contain. We considered three methods: a sequential method based on the analysis of increasingly complex molecules, an automated optimization method based on the bonds present in the molecules, and an automated optimization method based on the findings of the sequential method. The sequential method identified single N–O and N–N bonds together with NNO and ONO backbones as the largest error sources. On the other hand, the automated optimization method based on bonds deemed N–O, N=O, N–N, and O–H bonds as being problematic.

Comparison of the MAEs and MAXs among the three approaches indicates that the bond optimization method is the least accurate, while the optimization based on the errors detected sequentially is the best. This shows that (i) chemical intuition can be used to boost automated routines for error minimization, and (ii) the accuracy of the PBE functional to predict the thermochemistry of nitrogen compounds can be semiempirically enhanced, bypassing the need for more expensive levels of theory. Finally, we emphasize that the analysis shown here was made for PBE and nitrogen compounds but can be easily extended to other functionals and families of compounds.