Assessing Exchange-Correlation Functionals for Heterogeneous Catalysis of Nitrogen Species

Increasing interest in the sustainable synthesis of ammonia, nitrates, and urea has led to an increase in studies of catalytic conversion between nitrogen-containing compounds using heterogeneous catalysts. Density functional theory (DFT) is commonly employed to obtain molecular-scale insight into these reactions, but there have been relatively few assessments of the exchange-correlation functionals that are best suited for heterogeneous catalysis of nitrogen compounds. Here, we assess a range of functionals ranging from the generalized gradient approximation (GGA) to the random phase approximation (RPA) for the formation energies of gas-phase nitrogen species, the lattice constants of representative solids from several common classes of catalysts (metals, oxides, and metal–organic frameworks (MOFs)), and the adsorption energies of a range of nitrogen-containing intermediates on these materials. The results reveal that the choice of exchange-correlation functional and van der Waals correction can have a surprisingly large effect and that increasing the level of theory does not always improve the accuracy for nitrogen-containing compounds. This suggests that the selection of functionals should be carefully evaluated on the basis of the specific reaction and material being studied.


Introduction
Nitrogen is an essential component of all living things.Amino acids and nucleotides, the building blocks of protein and DNA, need nitrogen.Air is composed mostly of molecular dinitrogen, but the triple bond between two nitrogen atoms makes it inaccessible to use directly.Thankfully, certain bacteria can convert nitrogen gas molecules into useful compounds like ammonia through a process called nitrogen fixation.The converted nitrogen is used and transferred to other environments and organisms and then released as gas molecules at the end by the process called denitrification.Through these complicated processes, the nitrogen cycle is maintained and is critical for all living things.Human societies also require fixed nitrogen, particularly in the form of ammonia, for use as fertilizer.The Haber-Bosch process meets this demand and has significantly increased the world population.However, this has led to a significant artificial increase in fixed nitrogen, disrupting the natural cycle [1,2,3].Excess nitrogen fixation leads to the generation of greenhouse gases, acid rain, depletion of the ozone layer, and water pollution from nitrite ions [4,5,6].Thus, the management of the nitrogen cycle has been identified as a scientific "grand challenge" [7,8,9,10].Significant research has been devoted to overcoming this challenge in recent years, and heterogeneous catalysis, including electrocatalysis and photocatalysis, is a key strategy in efficiently converting nitrogen-containing compounds [11,12,13,14,15].
Computational studies are valuable for understanding and predicting heterogeneous catalytic materials.In the realm of nitrogen chemistry, numerous computational investigations have shed light on specific catalytic materials to synthesize and convert a range of nitrogen compounds [16,17,18,19,20,21,22,23,24,25,26,27].Ammonia synthesis is one of the most researched reactions in those studies.In computational studies, density functional theory (DFT) has been used to calculate the reaction rate on active catalysts [28,29], examine the uncertainty and reliability of computational results [20,30], and propose new catalyst design strategies [31,18,32] in ammonia synthesis.These successful applications of DFT to accelerate catalyst discovery and provide fundamental insight are not limited to ammonia synthesis.Computational techniques have been widely applied to study catalysts for other nitrogen-based reactions such as nitrate reduction [33,34], urea synthesis [35,36], and selective catalytic reduction [23,37,38,39,40,41,42], among others.
Discovering high-performance catalysts within the expansive realm of materials poses a formidable challenge.The volcano plot, rooted in the Sabatier principle [43], proves invaluable in identifying promising catalyst candidates by representing catalyst activity through adsorption energy descriptors [44].Taking ammonia synthesis as an illustration, a volcano plot can be constructed using the nitrogen adsorption energy as the descriptor [32,45].Researchers use this approach to find materials with optimal N adsorption energies to identify candidates with high catalytic activity and cost-effectiveness [46,47].As studies on ammonia synthesis have expanded, the scope of materials of interest (e.g., metal, metal oxide, nitride) and the range of reaction types (e.g., electrocatalysis and photocatalysis) have significantly increased.Exploring more intricate materials and reactions requires a mechanistic understanding at the atomic and electronic levels, a task achievable through computational analysis.An illustrative example is photocatalytic nitrogen fixation on TiO 2 (110) [48,49,50], where a combination of computational and experimental studies revealed the possibility of oxidative or carbon-assisted pathways, expanding the realm of possible mechanisms proposed from purely experimental studies.Intensive research is also underway on metal organic frameworks (MOFs) as promising materials for photocatalytic nitrogen fixation, with amine-functionalized MIL-125 serving as a notable example [51,52].Computational studies of nitrogen reactions on MOFs and other complex materials may yield valuable insights into the relevant mechanisms and strategies for improving catalyst performance.
Given the increasing importance of computational studies in nitrogen-containing compounds, ensuring the accuracy of computationally predicted properties is of utmost importance.DFT has become a prevalent tool in such studies and has been employed to compute the physicochemical properties of catalysts.While DFT provides computational efficiency compared to higher-level quantum calculations, it exhibits deviations in accuracy.This discrepancy arises from the inherent approximation within the exchange-correlation functional, a crucial element for describing electron interactions, leading to divergent computational results with different functionals."Jacob's ladder" of density functional approximation provides a general perspective on the cost-accuracy tradeoff of functionals [53].Those on the lower rungs (e.g., local density approximation, generalized gradient approximation (GGA)) are relatively fast (with cost typically scaling cubically with the number of atoms), but have larger errors on average.Functionals on the higher rungs (e.g., hybrids, random phase approximation (RPA)) are more accurate but can be drastically more computationally expansive due to quadratic or quintic scaling with the number of atoms, practically limiting their routine application to surface science systems.Nevertheless, the accuracy improvement only reflects a general trend, rather than a systematic guarantee, so the thoughtful selection of a functional tailored to a specific application is essential to achieve precise results with DFT.To address this, numerous studies have undertaken benchmarking of exchange-correlation functionals.Wellendorff et al. contributed significantly by releasing a benchmark database for adsorption energies on metal surfaces with six functionals [54] and conducting a comprehensive benchmark on various gas and solid phase properties, comparing the performance of BEEF-vdW with other functionals [55].Kepp conducted a specific benchmark for ammonia synthesis by modeling a single iron atom as a catalyst [30].Araujo et al. enhanced DFT predictions with corrections from small metal cluster calculations, improving accuracy in periodic systems, and benchmarked it on adsorption energies and activation barriers [56].Although there are various other benchmark studies [57,58], it is important to recognize significant limitations when extrapolating their results to reactions involving nitrogen-containing compounds.
One of the most critical constraints is the scope of materials and properties investigated in these benchmark studies.Metals dominate the benchmarked materials, with only a few exceptions.Other physicochemical properties have been benchmarked for metal oxides, but there are relatively few evaluations of chemisorption on oxides, particularly for nitrogen-containing species [59,60,61,62].Furthermore, a notable deficiency in most benchmarks is the absence of an adsorption energy comparison.Adsorption energy represents a fundamental property for the calculation of reaction descriptors, the illustration of reaction diagrams, and the prediction of reaction rates.However, it is notoriously challenging to benchmark due to the challenges in treating adsorption systems with high levels of theory or obtaining highly precise experimental measurements of chemisorption energies, particularly on complex materials [54].
In this work, we focus on selected nitrogen reactions, species, and solid catalyst materials relevant for nitrogen chemistry.We evaluated gas phase species relevant to a variety of nitrogen reactions (ammonia synthesis, nitrate reduction, urea synthesis, acetamide, nitromethane, and acetonitrile synthesis), and related adsorbed species (NO*, CN*, N*, N 2 H*, N 2 *, NH 2 *, and NH 3 *).These adsorbates are selected on the basis of their various types of covalent bonds (C−N, N−O, N−H, N−N) and their mechanistic importance in several previous studies of ammonia synthesis and nitrate reduction [63,64,28,49,40,33].For solid catalysts, we select two examples of metals (Pd, Cu), metal oxides (TiO 2 , MoO 3 ) and MOFs (MIL−125, OCUPUY).These materials are selected based on prior work on nitrate reduction and synthesis of ammonia or urea [65,40,66,49,48,67,68], existence of experimental adsorption energies [54], and presence of diverse bonding types including van der Waals layers (MoO 3 ) and open metal sites (OCUPUY).For functionals, we select several examples at GGA (PBE [69], RPBE [70]), GGA+vdW (BEEF-vdW [55], rev-vdW-DF2 [71]), meta-GGA (mGGA)(SCAN) [72], and hybrid (B3LYP [73], PBE0 [74], HSE06 [75]) levels of theory.We also utilize RPA optimized to gas-phase atomization energies to provide a ground truth for adsorption energies.These functionals are chosen on the basis of their common use in the literature and the diversity of physically derived and empirically fitted functionals.For all functionals without explicit vdW corrections, we also evaluate the influence of the empirical D3 parameters for dispersion forces.For gas-phase species, we compare to experimental formation energies, and for solid-state materials, we compare to experimental lattice constants.For adsorption energies, we compare to an optimized RPA result that, based on gas-phase performance, is expected to be close to chemical accuracy.Although the approach is far from exhaustive, it provides representative diversity in the electronic structures that are typical in nitrogen catalysis and the functionals commonly used to treat them, yielding the most systematic view of the role of exchangecorrelation functionals on nitrogen catalysis to date.
The results are roughly consistent with the conventional wisdom that functional choice leads to variations of ∼ 0.2-0.4eV in adsorption energies, with a standard deviation of 0.355 eV across all functionals, materials, and adsorbates (or 0.243 eV if the OCUPUY MOF with significant functional-dependent geometric differences is omitted).However, the errors compared to optimized RPA results are considerably larger (RMSE of 0.548 eV).The findings are also consistent with the expectation that the B3LYP hybrid functional has strong performance for gas phase formation energies (RMSE of 0.086 eV) and reaction energies (RMSE of 0.100 eV) [76].Interestingly, the B3LYP functional also yields excellent results for lattice constants, despite previous reports that hybrid functionals perform poorly for metals [77,78].We also find that the standard RPA functional yields results worse than B3LYP (RMSE of 0.110 eV), but that a few simple optimizations enable chemical accuracy (RMSE of 0.035 eV) on gas-phase formation energies with RPA correlation.Moreover, we find a surprising influence of dispersion corrections on adsorption energies, even with small molecules, leading to some counterintuitive conclusions.For example, as expected, the RPBE functional is one of the weaker binding functionals on average [54].However, when D3 corrections are included, it becomes the strongest binding functional, with D3 corrections that can exceed 1 eV, even for small adsorbates.We also find that, although the standard deviation between functionals is ∼ 0.2-0.4eV,specific deviations in the adsorption energies of individual species can vary by > 1 eV with different functionals.These deviations are not systematic across different nitrogen reactions, making it important to carefully select a functional based on the reaction of interest.

Details of electronic structure calculations
Density Functional Theory (DFT) calculations were performed using the Vienna ab initio Simulation Package (VASP) v6.1.2,employing projector-augmented wave (PAW) pseudopotentials (VASP v5.4.4 was used for BEEF-vdW calculations) [79,80,81,82].Dispersion corrections were incorporated in all calculations using the DFT-D3 method with the Becke-Johnson damping function [83,84], except for calculations that employ rev-vdW-DF2, BEEF-vdW or RPA, which inherently account for dispersion.The pre-defined damping function parameters for each exchange-correlation functional in VASP were used, with Grimme et al.'s parameters applied for HSE06 [85] (see Table S5).We omit +U corrections from all functionals due to ambiguities in the implementation of the method and the meaning of the U parameters across different implementations [86,87,88].
A k-point spacing of 0.4 Å−1 and a kinetic energy cutoff of 600 eV were used for bulk and slab models of Cu and Pd.For TiO 2 , MoO 3 , OCUPUY, and MIL−125, a k-points spacing of 0.5 Å−1 and a kinetic energy cutoff of 600 eV were used.The Γ point was included in all k-point grids.A kinetic energy cutoff of 600 eV was used to calculate the molecular formation energy and the molecular reaction energy.Convergence tests for the kinetic energy and k-points are provided in Table S1, Table S2, Table S3, and Table S4.Gaussian smearing of band occupancies with a smearing width of 0.05 eV was used.The convergence criterion of 10 −5 eV in total energy was applied to the self-consistent field (SCF) cycle.All calculations were spin-polarized.Nonspherical contributions related to the gradient of the density in the PAW spheres were included for hybrid functionals, SCAN, rev-vdW-DF2, and BEEF-vdW.The periodic and atomic cell relaxations continued until the forces on each atom were less than 0.03 eV/ Å.
Initially, each unit cell of bulk materials was relaxed to determine the optimized lattice constants and atomic positions.We used the initial structures of the two MOFs [89,90] from the CoRE MOF 2019 database [91].However, MIL-125 in CoREMOF 2019 database has deprotonated oxygen atoms, which is unphysical.We protonated the deprotonated MIL-125 as reported in prior papers [92,93].In the case of MOFs, the most stable spin state configuration was determined prior to lattice relaxation, with the tested spin states listed in Table S16.Using optimized cell and atomic positions, slab models for Cu (100 facet), Pd (111 facet), TiO 2 (110 facet), and MoO 3 (100 facet) were constructed in four layers, with the bottom two layers kept fixed in all subsequent calculations.To avoid self-interaction between repeated slabs, 10 Å of vacuum is added to the z-direction.For the slab models and the two MOFs, a gas molecule was introduced, and all atomic positions were optimized while the lattice constants was kept fixed.The initial placement of gas molecules was determined by manual placement on the basis of prior literature reports and intuition.The geometry optimization of adsorbates was performed solely for five non-hybrid functionals, and for three hybrid functionals and RPA functionals, single-point calculations were performed using PBE-optimized geometry.To obtain the energy of a molecule in DFT, a gas molecule was placed in the center of a cubic box with a 10 Å lattice and its geometry was optimized using a Γ point calculation.
All RPA calculations were performed using the low-scaling RPA algorithm [94] implemented in VASP v6.3.2.The GW pseudopotentials were used with 24 frequency points.For nonmetals, GW (H) or GW new (C, N, O) pseudopotentials were used, and for metals, sv GW pseudopotentials were used except for Cu, which was calculated using a GW pseudopotential with lower valency.The errors arising from pseudopotentials are expected to be < 0.05 eV on average (see Figure S5).For the calculation of DFT orbitals used in the evaluation of RPA correlation energy, the total energy was converged to below 10 −8 eV during the SCF cycle.
In a standard RPA calculation, the total energy is given by the equation, where E HF and E c,RP A refer to the Hartree-Fock energy and the RPA correlation energy, respectively, both evaluated non-self-consistently using DFT orbitals.The choice of DFT functionals leads to different orbitals and different RPA energies.RPA@PBE, RPA@PBE0, and RPA@PBEx50 refer to RPA calculations using different DFT orbitals obtained with PBE, PBE0, and a tuned hybrid functional, respectively.The tuned hybrid functional, PBEx50, is based on PBE0, but the amount of Fock exchange is set to 50 % as opposed to the 25 % used in PBE0.Previous reports have shown that including additional Fock exchange can improve the results of GW calculations [95].We also use an optimized version of RPA using the results of RPA with PBEx50 orbitals, where the RPA correlation energy is multiplied by a simple linear scaling factor of 1.17 (i.e.E optRP A = E HF + 1.17E c,RP A ) .The scaling factor for E c,RP A was determined to minimize errors in molecular atomization energies, as shown in Figure S1, based on the observation that the absolute contributions of E c,RP A to both the molecular atomization energy and the formation energy are underestimated in RPA calculations.It has previously been reported that the magnitude of E c,RP A can be incorrect [96], and applying the scaling factor to E c,RP A leads to highly accurate results for molecular energies.In the main text, we refer to this optimized version of RPA as optRPA, and it is always based on the PBEx50 orbitals.
For molecules, all RPA calculations were performed with a single Γ point, and molecules were placed inside cubic cells (with lattice lengths of 15 Å for E HF and 10 Å for E c,RP A ). Two different approaches were employed to calculate the gas formation/reaction energy and adsorption energy.To obtain the gas formation/reaction energy, the RPA energy was evaluated at the PBE optimized structures, and different kinetic energy cutoffs (700 eV for E HF and 600 eV for E c,RP A ) were used to obtain the two energy terms.For the adsorption energy, the relevant molecules were calculated with the same energy cutoffs as those used for the surfaces to maximize cancellation of error, and we verified that these lower cutoffs led to average numerical errors of <0.02 eV in molecular formation energies.
The energy of surface systems was determined at the PBE-optimized structures, with kinetic energy cutoffs of 600 eV for E HF and 400 eV for E c,RP A , based on the energy convergence test (see Figure S3).In the case of RPA calculations, varying formalisms, energy cutoffs, and smearing were used to achieve an expected numerical accuracy of 0.05 eV per system.Details of convergence tests for gas-phase, nonmetallic, and metallic systems are provided in the SI.

Evaluation
The evaluation metrics for the benchmarking comprised both gas-and solid-phase properties.The properties of the gas-phase molecules included molecular formation energy and molecular reaction energy.Molecular formation energy was computed using a least-squares regression approach (see Eq. (2) and Eq. ( 3)) to minimize systematic errors.Least-squares regression aligns calculated total energies of molecules (E i ) as closely as possible to experimental formation energy values (∆ f E i ) using atomic stoichiometries (n i,j ) and chemical potentials (µ j ) to minimize residual errors (ϵ i ).This approach ensures that all physically observable properties (i.e.energies of mass-balanced reactions) are consistent with the underlying method, but minimizes non-systematic cancellation of error.The resulting least-squares formation energies are denoted as ∆ LS f E i .Reaction energies were calculated using Eq. ( 4).All experimental energies of gas molecules were taken from the 0 K ATcT [97], and zero-point corrected using frequencies from the NIST Webbook [98] and CCCBDB [99].
For the solid-phase and interfacial properties, we computed the unit cell volume and gas adsorption energies.The volume of each DFT-optimized unit cell is compared to that of the experimental unit cell by normalizing the volume per atom.The experimental volume is preprocessed to directly compare it with the DFT-computed volume, by subtracting the thermal expansion volume and zero-point effect.Due to lack of experimental data, this thermal correction was not possible for MOFs, so the reported room-temperature lattice constants were used.The detailed methods and data used to correct the volume are in SI.The formation energy of the adsorbed species was calculated using Eq. ( 5), anchored to the reference species N 2 , H 2 , H 2 O, and CH 4 .∆ where E ref is a gas-phase reference based on N 2 , H 2 , H 2 O, and CH 4 .
Standard deviation evaluates the consistency between functionals in each crystal+gas system, except for RPAs.Total standard deviation (σ total ) of each material is computed by using Eq. ( 6).

Molecular Formation Energy
We selected ten nitrogen-containing species (CH ) commonly found in nitrogen-containing reactions (e.g., water splitting, CO 2 reduction).Figure 1 and Table 1 illustrate the deviation of calculated formation energy from the experimental formation energy.The optRPA functional has the smallest error in all metrics, and the RMSE of optRPA (0.035 eV) is within chemical accuracy (0.043 eV), and has a maximum error magnitude of 0.077 eV.Interestingly, the standard RPA (RPA@PBE) results are not significantly better than hybrid functionals (RMSE: 0.110 eV), and B3LYP-D3 exhibits the highest accuracy (RMSE: 0.086 eV) except for optRPA and RPA@PBE0 (RMSE: 0.054 eV), consistent with its well-known strong performance for gas-phase molecules [76].The other hybrid functionals, PBE0-D3 and HSE06-D3, yield very similar results to SCAN-D3, and actually exhibit slightly higher errors than SCAN-D3.At the GGA level of theory, BEEF-vdW yields the lowest errors, outperforming even SCAN-D3 and some hybrids, but notably, all 15 gas species were included in the training set for BEEF-vdW.The rev-vdW-DF2 and RPBE-D3 results are similar and slightly outperform PBE-D3, which yields the largest errors for gas-phase formation energies.

Molecular Reaction Energy
Reaction energies provide a practical metric for the accuracy of gas-phase properties, as they are directly relevant to chemical reactions.They differ from formation energies in that cancellation of error occurs implicitly in the reaction energy calculation and can thus yield different results from formation energy benchmarks.
Here, we select six common reactions in nitrogen chemistry.The difference between the calculated and experimental reaction energies for these reactions is depicted in Figure 2, and the corresponding errors for each functional are presented in Table 1.Trends are broadly consistent with the formation energy results, with optRPA and B3LYP-D3 emerging as the most accurate overall.Interestingly, B3LYP-D3 (RMSE: 0.100 eV) outperforms RPA@PBE0 (RMSE: 0.128 eV) in the reaction energy, mainly due to large error of RPA@PBE0 for NO reduction to NH 3 .optRPA exhibits errors slightly higher than chemical accuracy on average (RMSE: Notably, these conclusions may differ if other reactions are selected, but the focus here is on evaluating common nitrogen chemistry reactions.One advantage of the BEEF-vdW functional is its ability to provide error estimates.These error estimates reflect the sensitivity to the GGA exchange enhancement factor and are generally expected to reflect variations that arise from selecting different GGA functionals.However, it is clear that this is not always the case, with RPBE-D3 and PBE-D3 often falling well outside the error bars of BEEF-vdW.This may be due to effects of different geometries between the functionals, the fact that BEEF-vdW error estimates are not based on self-consistent densities, or the inclusion of D3 corrections in other GGA functionals.Other classes of functionals, particularly hybrids and RPA functionals, are often even further outside the error estimates of BEEF-vdW [103], indicating that BEEF-vdW error bars should be interpreted with caution.Nevertheless, there is a general correlation between the BEEF-vdW error estimate and the standard deviation of all functionals for a given reaction, suggesting that, although the error estimates are not quantitatively accurate, they reliably capture the relative sensitivity to functional choice.
Another interesting observation from the reaction energies is the high accuracy of reaction energies for NO synthesis from N 2 and O 2 for all functionals.This is remarkable, given the fact that O 2 has a triplet spin state and that NO contains an unpaired electron.These complex electronic structures are generally Table 1: Error of the calculated formation energies and reaction energies compared to the experimental formation and reaction energies.Unit in eV.difficult to treat with DFT, and O 2 is commonly "corrected" by 0.4 -0.5 eV to ensure the correct energetics of the water splitting reaction [104].The high accuracy of NO synthesis with all functionals, along with the relatively large errors for other reactions involving NO and O 2 , suggests that there is reliable error cancellation between O 2 and NO, and reveal that care should be taken if any O 2 correction is applied in cases where NO synthesis also occurs.

Lattice Constants
Unit cell volume is selected as a metric for assessing solid-state performance of functionals, because experimental values are generally available for all materials classes of interest.A comparison between the calculated and experimental unit cell volume per atom and errors is shown in Figure 3 and Table 2.The OCUPUY MOF exhibits considerable variance because of the potential for multiple geometries of the water molecules in the structure.This leads to inconsistent optimization of the unit cell and atom positions, with SCAN-D3 and PBE0-D3 predicting different geometries of water than other functionals.These effects are difficult to deconvolute, as discussed in the Discussion, so OCUPUY was excluded from the computation of max error, MAE, and RMSE to provide a more fair comparison of performance.Across the five materials, B3LYP-D3 shows the smallest RMSE, even for metallic systems, contrary to some previous reports on the applicability of B3LYP to solid and metallic systems [77,78].The HSE06-D3, PBE0-D3, and SCAN-D3 functionals also have comparable accuracy to B3LYP-D3.Of the GGA functionals, rev-vdW-DF2 and PBE-D3 demonstrate relatively high accuracy, while BEEF-vdW and RPBE-D3 have substantially larger errors.The low accuracy of BEEF-vdW may be due to the fact that oxides and MOFs were not included in the training set, although the error on metals is also higher than that of most other functionals.The poor performance of RPBE-D3 is attributed to the D3 corrections, as discussed further in the Discussion.We note that the very small number of solid-state materials included here and the fact that only cell volume is evaluated makes it difficult to draw strong conclusions about the specific ranking of functionals for solid-state properties.However, the diversity of materials classes evaluated is expected to yield general qualitative insight into the suitability of these functionals for the treatment of solids.

Adsorption Energies
Adsorption energies for several nitrogen-containing intermediates (N*, N 2 *, CN*, NO*, N 2 H*, NH 2 *, NH 3 *) calculated at different levels of theory on all six materials are illustrated in Figure 4. Given the absence of experimental adsorption energy data for most systems, we utilize the optRPA result as the "ground truth" on the basis of its excellent performance for gas phase energies (see Table S9 for the calculated adsorption energies from optRPA together with other RPA variants).The commonly used RPA@PBE approach differs from RPA@PBE0 by an RMSE of 0.217 eV, and from optRPA by an RMSE of 0.297 eV for systems where all RPA results are available (all Cu systems and NH 2 , N 2 H, and CN on TiO 2 ).The errors of all other calculated adsorption energies compared to optRPA are listed in Table 3 and Table 4. Calculated adsorption energies from all four RPAs (optRPA, RPA@PBEx50, RPA@PBE0, RPA@PBE) are listed in Table S9.For systems with known experimental values (Cu(100)+NH 3 *, Pd(111)+NO*), the experimental adsorption energy is marked as a yellow star, and is very close to the optRPA result.As discussed in the Methods, all hybrid and RPA calculations are obtained through single-point calculations using PBE-D3 optimized geometry.The numbers in this section represent a practical comparison of adsorption energies that would be obtained if standard practices in the literature were followed; a more detailed discussion of the influence of geometries and spin states is provided in the Discussion.
In general, the deviations of adsorption energies between functionals for metals and oxides exhibit standard deviations ∼0.2 eV, consistent with the commonly accepted exchange-correlation error of 0.2 -0.3 eV.The standard deviations of all functionals for metals and metal oxides are 0.251 eV and 0.274 eV respectively.Interestingly, if optRPA is assumed to be the ground truth, metals have a higher RMSE (0.654 eV) compared to metal oxides (0.365 eV).All functionals systematically overestimate the adsorption strength in metals with an average mean signed error (MSE) of −0.522 eV.There is also an overestimation trend in metal oxides, but the amount of overestimation is much smaller (MSE: −0.134 eV).On the other hand, RPA@PBE0 tends to bind weakly, with MSEs ranging from −0.025 eV on metals to 0.128 eV on MOFs, and exhibits RMSEs of 0.143 eV for metals, 0.157 eV for oxides, and 0.233 eV for MOFs.Caveats regarding optRPA adsorption energies are discussed in the Discussion.For metals, the BEEF-vdW functional has the lowest adsorption energy error (RMSE: 0.447 eV) compared to optRPA, which is consistent with the fact that it was trained on metal adsorption energies, although very few of the systems evaluated here were included in the training set [55].Similar to the case of gas-phase energies, the BEEF-vdW error bars generally do not encapsulate other functionals, although their size loosely correlates with the spread of all functionals with a Pearson correlation coefficient of 0.536 (see Figure S12).In metals, hybrid functionals have performance similar to BEEF-vdW and each other, with PBE0-D3 and B3LYP-D3 being slightly more accurate than HSE06-D3.On the other hand, all other GGA functionals and SCAN strongly over-bind, with rev-vdW-DF2 having the best performance, and PBE-D3 being slightly worse.Surprisingly, RPBE-D3 shows the worst performance and significantly overbinds in all cases, which is true even in metal oxides except for a few instances on MoO 3 .This is contrary to the expectation that RPBE binds relatively weakly and is due to D3 corrections, as discussed further in the Discussion.Interestingly, the tendency to overbind is stronger on metals than on metal oxides for all functionals, and all functionals are more accurate for oxide adsorption energies than metals, despite the fact that no +U corrections are applied.However, for metal oxides the relative performance of functionals shifts, and hybrids are more accurate than BEEF-vdW.In TiO 2 , B3LYP-D3 exhibits the lowest error among the eight functionals (RMSE of 0.194 eV).The PBE-D3, PBE0-D3 and rev-vdW-DF2 functionals exhibit relatively little systematic error on oxides (MSE < 0.1 eV), and rev-vdW-DF2 performs particularly well (RMSE of 0.211 eV) for the MoO 3 system with vdW layers.RMSE for each material is listed in Table S11.
In the case of MOFs, the errors and deviations are large because of OCUPUY, with a total standard deviation between all functionals of 0.489 eV, and a RMSE relative to optRPA of 0.583 eV across all functionals.The situation is also quite different between the two MOFs chosen in this study.In the case of OCUPUY, there is relatively little systematic error compared to optRPA (MSE = −0.041eV), but there is a very large deviation between the functionals (standard deviation of 0.678 eV).On the other hand, for MIL-125 there is a relatively larger systematic error compared to optRPA (MSE = −0.183eV), but a smaller deviation between functionals (standard deviation of 0.140 eV).This reflects a key difference in the two MOFs, since adsorption occurs at a spin-polarized open metal site in OCUPUY, but occurs in the pores of MIL-125.Furthermore, as discussed in the results of lattice constants, OCUPUY contains several weakly bound water molecules that adopt different geometries for different functionals.The effects of geometry and magnetism are discussed in more detail in the Discussion, but here we focus on the deviations that would be observed in a typical high-throughput study with standardized settings across materials.In the case of MIL-125, the adsorption energies are much more consistent between functionals, with a small standard deviation of 0.140 eV between functionals, and a maximum error of ∼0.5 eV when compared to optRPA.Among non-RPA functionals, B3LYP-D3 has the smallest error (RMSE of 0.169 eV).The hybrid functionals and the non-hybrid functionals have similar errors in MIL-125 while the hybrid functionals are much more accurate in OCUPUY.This is consistent with previous reports that hybrid functionals required for open metal sites [88,105].These findings indicate that chemisorption in MOFs with open metal sites may exhibit errors much larger than those of metal or metal oxide systems.We found that the source of the large error includes not only differences between functionals but also variations in geometry, though small amounts of magnetism also contribute.The effects of different geometry and magnetism are discussed in the following Discussion section.For Pd, the adsorption energy from RPA@PBE is set to be the ground truth because optRPA calculations require significantly higher computational cost due to the use of the densest k-point grids.The impact of the D3 correction on calculated adsorption energies was found to be substantial, showing notable variations between different functionals and systems.To examine the contribution of the D3 correction, we compared properties with and without the correction.The latter were obtained by subtracting the D3 correction from the D3 corrected results.As expected, molecular formation energies and molecular reaction energies are relatively insensitive to the D3 correction (see Table 1, Table S7, Figure S7, and Figure S8), although the D3 correction has a noticeable effect on the gas-phase properties calculated by RPBE-D3.This is because the parameters of RPBE (s 8 , a 1 , and a 2 ) used in Eq. ( 8) for the D3 correction make the calculated dispersion energy large.Interestingly, almost all functionals have slightly increased RMSE and MAE for gas phase properties when the D3 correction is added.However, the difference is ∼0.02 eV in RMSE, which is comparable to the numerical accuracy of the calculations.
In contrast, the optimized lattice constants are significantly impacted by dispersion correction.To test the effect of D3 correction on unit cell volume, we performed cell relaxations without the D3 correction on metals and metal oxides using PBE (see Table S8).Due to the absence of dispersion interactions, cell relaxation without D3 correction results in expanded volume compared to Figure 3.The amount of expansion is significant, so the discrepancy from the experimental volume becomes larger.In particular, the lattice constant of MoO 3 , a material with van der Waals layers, changes drastically (∆V = 2.215 Å3 /atom) without D3 correction.The change was even greater with RPBE (∆V = 3.288 Å3 /atom) due to the larger magnitude of the parameters used in Eq. ( 8).The conclusion is that the D3 corrections generally improve the volumes of solids.
Despite relatively small changes in gas phase energies and improvements in solid-state lattice constants, the effects of D3 corrections on adsorption energies are significant, highly varied, and difficult to assess.Surprisingly, the change in adsorption energy after adding the D3 correction can exceed 1 eV depending on the systems and functionals.Because the D3 correction adds attractive dispersion interactions, adsorption becomes stronger when the D3 correction is included in all cases except for NH 2 adsorption on MoO 3 using RPBE-D3.This exception occurs due to reconstruction of atoms in MoO 3 after the adsorption of NH 2 .Figure 5 shows the distribution of D3 corrections as a function of the adsorption energies, with box plots showing deviation by functional and adsorbate.Interestingly, there is only a very weak correlation between the magnitude of the D3 correction and the magnitude of adsorption or the size of the molecule, revealing that contrary to conventional wisdom, the D3 correction can have a significant influence even in cases of strong chemisorption of atoms or small molecules.
There is also a significant variation in the magnitude of D3 corrections between functionals.In the case of RPBE, the magnitude and variance of the D3 correction was much larger than other functionals, with an average contribution of −0.743 eV, causing RPBE-D3 to be one of the strongest binding functionals.However, if the D3 correction is removed, then RPBE becomes one of the weakest binding functionals (Figure S10), in agreement with prior work [54].Based on available experimental adsorption energies (Table S13) and optRPA results, it appears that the D3 correction overcorrects the adsorption energies for RPBE.B3LYP has the second largest D3 correction, −0.470 eV on average, which yields more accurate results for the limited number of experimental adsorption energies.The D3 corrections for PBE, PBE0 and HSE06 are smaller in magnitude and similar to each other (∼−0.29 eV on average), while SCAN has the smallest correction (−0.103 eV on average).The corrections for these functionals are more consistent with the typical magnitude of physisorption, and give adsorption energies similar to rev-vdW-DF2, so they are likely more accurately capturing dispersion in the adsorption systems.
Notably, in VASP the D3 correction also influences both energies and forces and can significantly affect the geometry of the system, especially when dispersion interactions are crucial to maintain the crystal structure.To deconvolute the effects of geometry and D3 dispersion energy, it is necessary to compare the adsorption energy with the interaction energy, which is the energy difference with all geometric coordinates fixed [106].We compare these two quantities for three representative cases: (i) CN+MoO 3 with strong chemisorption and dispersion interactions in the solid, (ii) CN+Pd (111) with strong chemisorption and no significant dispersion forces in the solid and (iii) N 2 +Pd(111) with dispersion forces dominating the adsorption of the weakly bound adsorbate (see Table S14).In the case of MoO 3 +CN, the geometry changes played a significant role.The D3 correction in the interaction energy (−0.57eV) is much smaller than in the adsorption energy (−1.46 eV), indicating that the dispersion-induced reconstruction of MoO 3 contributes significantly to the D3 correction.However, in the cases of N 2 and CN adsorption on Pd (111), the contribution of D3 is independent of geometry.This indicates that the dominant effect in most cases is the energetic contribution of the D3 correction, rather than the geometry changes.

Influence of Geometry and Spin States
Adsorption energies are the critical parameter in computational studies of heterogeneous catalysis, and are hence the primary quantity evaluated here.However, while properties such as unit cell volume and gas state energies are straightforward to compute with a single DFT simulation, calculating adsorption energies requires multiple steps and presents a more intricate challenge in determining the accuracy and precision of each functional.An adsorption energy calculation involves lattice relaxation to find the optimal unit cell of the solid, computing the energy of the gas phase molecule (or reference), computing the energy and geometry of the adsorbent without the adsorbate present, and finally calculating the energy of the adsorbateadsorbent system.Each of these steps involves an optimization of both the atomic positions and magnetic moments (if the system is magnetic), both of which may differ between functionals and/or between the adsorbent systems with and without the adsorbate present.Differences in geometry or magnetic moments between any of the subsystems involved will be included in the adsorption energy.These differences in geometry and magnetic moments may reflect real physical differences in the system but may also lead to artifacts where not all systems are in their global minimum state.This issue is exacerbated in comparison of functionals, where different functionals may have significantly different geometric and/or magnetic ground states.The adsorption energies presented in the Results section encompass all of these effects and are meant to provide practical insight into the deviations in adsorption energies between functionals in highthroughput scenarios where standard input settings are used.In this section, we briefly outline scenarios where the indirect geometric and/or magnetic effects may be as or more significant than the direct impact of the energy functional.
Different exchange-correlation functionals always have slight discrepancies in geometric configurations, but in some cases they lead to entirely different geometries that are not equivalent.For example, NO adsorption on MoO 3 exhibits two distinct local minima, chemisorption and physisorption, where specific GGA functionals (PBE-D3, RPBE-D3) capture both states while others converge only to chemisorption regardless of the initial guess used (see Table S15).Additionally, the relaxation of OCUPUY's cell reveals varied geometries for water molecules binding to vanadium atoms depending on the functional used, significantly impacting the optimized volume of OCUPUY.Retaining solvent molecules affects the chemical environment of the open metal site, influencing electronic properties and applications such as gas adsorption.Often computational studies are conducted on the clean MOFs by removing all solvent molecules perfectly for consistency and simplification, but this can lead to atomic geometries of open metal sites that are not achievable in real experiments.The process also commonly leads to artifacts such as deprotonated linkers that can drastically influence results.Notably, these artifacts are present even for well-established databases such as CoRE MOF 2019 and common MOF materials such as MIL-125.Figure S11 shows the drastic (> 3 eV) differences in adsorption energies that arise when the deprotonated structure from CoRE MOF 2019 is used directly, compared to the properly protonated structure.Furthermore, the original OCUPUY paper reported water crystallization in the framework [90], so solvent molecules were retained in OCUPUY in this work.Lattice optimization of OCUPUY yields two distinct local minima (Figure S9), primarily due to variations in the orientation of the water molecule bound to vanadium open metal sites, leading to significant variance among functional calculations and introducing errors not solely attributed to the accuracy of the exchange-correlation functional.Despite significant effort, it was not practically feasible to force OCUPUY to conform to a consistent geometry between functionals during the lattice optimization, highlighting that these geometric effects are often unavoidable in practice.
In addition, the transition-metal atom at open metal sites in MOFs often have unpaired electrons, resulting in multiple magnetic states (high-spin or low-spin).Within a MOF unit cell, each metal atom possessing unpaired electrons can adopt different magnetic configurations, either up or down, leading to a combinatorial problem to determine optimal spin configuration.In high-throughput settings, this optimization is typically performed for the empty MOF structure [106], but magnetic configurations can change during gas molecule adsorption.For example, in OCUPUY, we verified that the optimal spin state configuration comprises all Figure 6: Error of the calculated adsorption energies in OCUPUY from the adsorption energy of optRPA calculation, with the (left) fixed PBE-optimzed geometry and PBE-optimized magnetic moments, (center) the fixed PBE-optimized geometry but relaxed magnetic moments, and (right) all relaxed systems (same with the OCUPUY results in Figure 4).PBE Geometry refers to geometry from PBE-D3 optimization.up spin states for the vanadium atoms (Table S16), yet this optimal state is not always retained during gas adsorption processes (Table S17).In numerous computational studies, the focus has been on determining the optimal spin state of the MOF adsorbent [107,108].However, the majority of high throughput studies have not accounted for potential changes in the optimal spin state when gas molecules are adsorbed.Therefore, Figure 4 represents the variations that are likely present in results generated with high-throughput computational workflows.Similarly to the case of geometry differences, these changes in magnetic moment may differ between functionals, further convoluting the effect of the exchange-correlation functional.Notably, Hegde et al. reported that magnetization is one of the most difficult properties to reproduce across all classes of solid-state materials, even when working with the same input parameters as the database that reports these properties [109].This indicates that the adsorption energies for magnetic systems will generally exhibit much larger variations between different functionals and different studies.
As a case study in the influence of geometry and magnetism, we deconvolute these effects in the case of adsorption on OCUPUY, since it shows the largest deviation in adsorption energies and contains significant differences in geometry and magnetic states between functionals.As a most constrained system, we conducted single point calculations on PBE-optimized geometry and PBE-optimized magnetism (using the PBE-opimized wavefunction as an initial guess) with different functionals (PBE-optimized refers to PBE-D3optimized in this section).In addition, to solely evaluate the effect of different magnetization in adsorption energy, single point calculations on PBE-optimized geometries without using the PBE-optimized magnetism, but setting +1µ B to all atoms (the default setting of VASP when running a spin-polarized calculation) are conducted.The results are shown in Figure 6 along with the original adsorption energy result from Figure 4.In addition, the difference in adsorption energy, relative to the most constrained "PBE Geometry & PBE Magnetism" system is provided in Table 5.When geometry and magnetic moment are fixed, the variation between non-hybrid functionals is similar to what is seen in metals and metal oxides, with a standard deviation of 0.326 eV.The effect of magnetism alone is relatively minor in OCUPUY, with the standard deviation between functionals increasing to 0.350 eV when magnetic moments are relaxed, although in certain cases such as N 2 H the influence can be much larger for some functionals.Relaxing the geometric constraints is the dominant effect in the variation between functionals, leading to a total standard deviation of 0.662 eV between non-hybrid functionals.In the case of N 2 H adsorption in BEEF-vdW, the difference between original adsorption energy and the one from fixed PBE geometry and PBE magnetism reaches nearly ∼1.5 eV.SCAN-D3 shows an even larger discrepancy when relaxing geometry, which arises due to the different unit cell geometry caused by the different orientation of the water molecule.In summary, systems like MOFs with many geometric degrees of freedom and multiple possible magnetic states are likely to exhibit very large

Assessment of Functionals for Nitrogen Catalysis
While many results from this work hint at general trends in accuracy of exchange-correlation functionals, the primary scope is to evaluate functional choice for nitrogen reactions.Unfortunately, there is no clear "best choice" of functional, and the selection will need to be based on a tradeoff of cost, accuracy, specificity, and generality.If a general-purpose functional is needed that works across many nitrogen reactions and catalyst types, then B3LYP-D3 is a strong choice if hybrid functionals are not too expensive.It demonstrates low error across all gas-phase reactions (unlike HSE06-D3 and PBE0-D3 which have errors ∼0.4 eV for NH 3 reactions), good solid-state lattice constants, and competitive adsorption energies compared to RPA.If studying ammonia or NO synthesis from N 2 on metals, the BEEF-vdW functional is an excellent choice that gives strong gas-phase reaction energies and adsorption energies, but if other reactions such as NH 3 oxidation are of interest, then gas-phase errors become problematic.The SCAN-D3 and rev-vdW-DF2 functionals are also good general-purpose choices with comparable accuracy across most metrics.Despite its common use, the PBE-D3 functional exhibits relatively high errors for properties other than adsorption on oxide materials, and the RPBE-D3 functional shows the largest overall errors, although this is attributed to D3 parameters that are not tuned for adsorption.There are a few specific nitrogen adsorption systems that also warrant discussion.Amine-functionalization of MIL-125 is common in CO 2 capture [110,111,112] and has been reported to be an active photocatalyst for nitrogen fixation [68].Our adsorption energy results show that B3LYP-D3 has the smallest error, while the other two hybrid functionals (PBE0-D3 and HSE06-D3) and vdW functionals (rev-vdW-DF2 and BEEF-vdW) also exhibit small error.Regarding the accuracy of gas adsorption energy in MIL-125, B3LYP-D3 appear to be the best choices among hybrid functionals because accuracy of HSE06-D3 and PBE0-D3 in the ammonia synthesis reaction is low, so B3LYP-D3 may be a better choice since it exhibits reasonable accuracy in both adsorption energy and reaction energies.In addition, rev-vdW-DF2 and BEEF-vdW can be accurate because of the small error in adsorption energy (RMSE of 0.261 eV and 0.238 eV respectively) and ammonia synthesis reaction (error of −0.055 eV and 0.109 eV respectively).Another important case is TiO 2 , which has been extensively studied for photocatalytic ammonia synthesis [48,113,114,115,116,117].In most cases, the BEEF-vdW adsorption energies are similarly accurate to hybrid functionals, although there is a substantial deviation for CN adsorption, where the error of BEEF-vdW and other GGA functionals is ∼1 eV.This suggests that hybrid calculations may be necessary for evaluating carbon-assisted pathways [49,50], and that the HSE06-D3 or PBE0-D3 adsorption energies of NH x species can be considered accurate despite the large error in the gas phase ammonia synthesis reaction.
There are also several notable observations about NO.As discussed in 3.2, the reaction energy for NO synthesis from N 2 and O 2 is remarkably insensitive to functional choice, and is at or close to chemical accuracy for any selected functional.This is surprising, given the complexity of the electronic structure of the O 2 and NO molecule.However, adsorption energies of NO exhibit much larger errors and variations, so functional choice is still important for the NO synthesis reaction.NO adsorption energy errors on Pd and Cu are ∼ 1 eV on average, with a significant spread (std dev = 0.366 eV) between functionals.The situation is somewhat better for other materials, but large deviations between functionals are observed.This, along with the well-known challenges in treating O 2 [104], indicates that the accuracy of gas-phase NO synthesis is not due to an accurate treatment of NO or O 2 by the functionals, but rather a surprisingly consistent cancellation of error in all cases.The electronic structure origins of this phenomenon are beyond the scope of this study, but may yield useful insight into nitrogen oxidation or exchange correlation design in future work.

Limitations of Dataset and Methods
This study provides deep, but narrow, insight into the choice of exchange correlation functionals for nitrogen chemistry in heterogeneous catalysis.The dataset is one of the most diverse available for comparison of adsorption energies at the RPA level of theory, but it is very sparse compared to the possibilities of adsorbates and catalyst materials, and it is intentionally biased toward nitrogen chemistry.It is difficult to draw general conclusions from the results of the six materials and seven adsorbates studied here, especially given the chemical diversity of the materials involved.It is possible that the materials selected are not actually representative of metals, metal oxides, or MOFs, and hence additional tests should be performed on any given system of interest before drawing strong conclusions about the accuracy of a particular exchange correlation functional.However, some clear trends emerged from the small and diverse dataset studied here, and even the results from the handful of selected materials provide insight into some significant pitfalls and unexpected outcomes in the performance of exchange correlation functionals in nitrogen chemistry.
The lack of a well-established ground truth for the adsorption energies is also a limitation of this work.Here, we selected the optRPA level of theory as a ground truth based on the level of physics included and its strong performance in gas-phase energetics.However, we note that, to our knowledge, the particular strategy of optimizing the correlation contribution has not been previously applied, so while the gas-phase results are promising, the approach is not well-tested outside the scope of the gas-phase reactions and atomization energies studied here.The effect of the optimization of the RPA correlation can be estimated by comparing the RPA@PBE0 and optRPA results, and is relatively small in most cases, so we expect that the transferability of optRPA from molecules to solids is generally reliable.Another possible source of bias is the selection of wavefunctions used in optRPA, which are closely related to the PBE0 wavefunctions.This may bias the accuracy of hybrid functionals toward PBE0, although the modification of the exact exchange contribution from 25% to 50% is significant, and the results of the atomization energy benchmarks suggest that the amount of exact exchange included is more significant than the functional choice (see S1). Fully self-consistent RPA correlation is required to evaluate this effect, but is beyond the scope of this study.
An additional limitation is that the RPA calculations are far more involved than standard DFT simulations.This human cost, along with computational cost, is illustrated in Figure S14.RPA calculations require testing convergence of parameters like electronic smearing that are less important in standard DFT.The RPA correlation is calculated separately and converges differently from the exact exchange, necessitating a careful combination of calculations performed with different settings.Details of all RPA calculations are provided in the SI, and we expect that they are converged with a numerical error of ∼0.05 eV on adsorption energies.Notably, this required extreme values of k-point densities and electronic smearing, especially in the case of metals, leading to massive computational costs.Failure to properly converge any of the numerous parameters involved in RPA can easily lead to numerical errors that exceed 0.5 eV (see Figure S4), so it is possible that RPA numbers obtained from the literature (or even those in this work) may have uncontrolled numerical errors.However, strong agreement with gas-phase properties and extensive convergence testing provides confidence that the numbers included here are sufficiently accurate to be considered a ground truth.Notably, we find that the computational cost of RPA single point calculations is not prohibitive for typical adsorption systems, but memory requirements are higher, and significantly more convergence testing and manual analysis is required to achieve reliable results, as illustrated in Figure S14.It is also worth noting that RPA using wavefuctions from global hybrid functionals is numerically more stable compared to the RPA@PBE approach.This indicates that optimized RPA can be selected when higher accuracy is needed than hybrid functionals, but extra caution in choosing parameters for the calculation to obtain reliable and converged results.

Conclusion
In the present work, we benchmarked multiple functionals from vairous categories (GGA, meta-GGA, hybrid, RPA) for nitrogen chemistry, using gas-and solid-phase properties.The results indicate that gas-phase errors for nitrogen reactions can be very large (>0.4 eV) even for some hybrid functionals, but that with a few minor modifications it is possible to achieve near-chemical accuracy in gas phase properties at the RPA level of theory.We utilized this optimized version of RPA to provide highly accurate benchmark adsorption energies of nitrogen species across a diverse range of solid materials: metals, metal oxides, and MOFs.The results reveal that, although standard deviations of adsorption energies between functionals fall roughly in the 0.05 -0.30 eV range, the RMSEs when compared to accurate benchmarks can be considerably larger (>0.6 eV).Surprisingly, the adsorption energy errors were lowest on average for the oxide materials, despite the presence of highly localized electrons and van der Waals layers.The results generally reveal tradeoffs between the functionals studied, with B3LYP-D3 demonstrating a good compromise between gas, solid, and adsorption properties across all systems studied.At lower levels of theory, the BEEF-vdW functional showed excellent performance for metals, but had large errors for gas-phase reactions or materials outside its training domain.The rev-vdW-DF2 functional also showed strong performance for gas, solid and adsorption properties, though, like most nonhybrid functionals, it exhibits strong overbinding on metals.Overall, the small number of materials studied here makes it difficult to draw strong conclusions but provides some insights into which functionals to more carefully evaluate for a given material and reaction in nitrogen catalysis.
The work also revealed several interesting insights that should be generally considered when evaluating adsorption energies with DFT.Empirical D3 corrections can vary wildly between different functionals and in several cases lead to dispersion contributions that exceed 1 eV, even for small and strongly chemisorbed molecules.This is particularly true for RPBE, where the D3 correction takes it from one of the weakest binding to one of the strongest binding functionals.Additionally, we find large deviations in the adsorption energies of MOFs when open metal sites are present.These deviations can be largely attributed to changes in geometry, indicating that the flexibility of MOFs makes calculated adsorption energies more sensitive to functional choice.The presence of different spin states in open metal sites also contributes to these deviations, but is less significant than the geometry in the case studied here.However, in the case of MIL-125, small deviations in adsorption energies (standard deviation of 0.140 eV) were commonly observed between functionals, and RMSE of 0.270 eV was observed compared to the optimized RPA results.Relatively little effort has been applied to using RPA or other double-hybrids on MOFs due to their large system sizes, but these results suggest that further benchmarking and analysis is needed to reliably determine when standard DFT functionals can accurately predict the adsorption of nitrogen species in MOFs with open metal sites.
The present work provides practical insight into the variations expected between different functionals when studying nitrogen catalysis, and sheds light on the effects of multiple artifacts that researchers will encounter by following the standard workflow of adsorption energy calculations in studies of heterogeneous catalysts.Even with the small number of materials and adsorbates studied, it is clear that all available functionals have some weaknesses, and that care should be taken in selecting a functional if quantitative accuracy is desired.We hope that this work serves as a useful guide for computational researchers in nitrogen catalysis who are faced with the problem of functional selection.approach.The error of optRPA@PBEx50 is comparable to the numerical error in CCSD(T), as can be seen in the error of CCSD(T)-F12, a CCSD(T) method with correction, which improves the basis set dependency.[126] The same scaling approach would not work for RPA@PBE as we have seen examples where the signs of E HF and E c,RP A of molecules were not correct.In the assessment of the sensitivity of the E c,RP A scaling factor using k-fold cross-validation, we observed consistent values for the scaling factor in different fold configurations.For the optRPA@PBEx50 approach, considering k=3, 5, and 10 folds, the mean (median, standard deviation) of the optimized scaling factors were found to be 1.168 (1.168, 0.0002), 1.167 (1.168, 0.0010), and 1.167 (1.168, 0.0007), respectively.The scaling factor optimized based on the formation energy was identified as 1.174.When employing B3LYP-optimized structures and atomization energies, and the scaling factors were determined to be 1.176 for optRPA@PBEx50 and 1.248 for optRPA@PBEx75.B3LYP-optimized structures were used to obtain the energies for Figure S1 and Figure S2, but deviations when PBE geometries are used were observed to be minimal in all cases tested.
RPA Energy Convergence Detailed convergence test results are provided for RPA calculations across various systems.RPA calculations are computationally intensive, especially in terms of memory requirements.Therefore, a lower kinetic energy cutoff (ENCUT) of 400 eV was chosen for calculating ∆ r E c,ads as defined below, which represents the E c,RP A contribution to the adsorption energy ∆ r E ads , after a convergence test (Figure S3).
∆ r E ads = E solid+gas − E solid − E gas (S3) Among the GW pseudopotentials used, the oxygen (O GW new) and nitrogen (N GW new) pseudopotentials have the highest default ENCUTs of 434.4 and 420.9 eV.Consequently, the observation that, for MoO 3 systems containing N-containing adsorbates, the differences in ∆ r E ads between ENCUTs of 400 and 500 eV are less than 0.02 eV suggests that 400 eV is a reasonable choice (Figure S3(a)).Although differences of about 0.05 eV are observed for OCUPUY, the absolute ∆ r E c,ads exceeds 5 eV, corresponding to roughly a 1% error (Figure S3(b)).
In the RPA@PBE approach to metallic systems, the non-finite-temperature formalism may break down; thus, the finite-temperature formalism should be used.Additionally, the calculations on metallic systems are sensitive to σ, as the metallic slabs, with and without adsorbates, respond differently to variations in σ (Figure S4).The different responses can be somewhat alleviated by subtracting half of the electronic entropy term from E HF , similar to how the zero-temperature energy is obtained in DFT calculation.Furthermore, a correction to E HF related to partial occupancies [127] was not applied, as it worsens the energy convergence in terms of σ (Figure S4(a)).According to the VASP forum, the correction cannot be rigorously derived within the finite-temperature formalism.Therefore, the correction to E HF was applied only to calculations using the non-finite-temperature formalism.RPA calculations on metallic systems generally require much denser k-point grids than standard DFT.For the metallic slabs of Cu and Pd, k-point grids of 7×7×7 and 12×12×12 were used, respectively (See Figure S4(b) and Figure S5(b) for energy convergence test).Another parameter worth discussing is the choice of GW pseudopotentials.Specifically, for Cu, a GW pseudopotential with fewer valence electrons was utilized after convergence test (Figure S5(a)).The reason behind the choice is the high computational cost associated with generating PBEx wavefunctions, rather than the actual RPA calculation step to get E c,RP A .
RPA@PBE calculations on TiO 2 reveal a sensitivity to σ of non-metallic systems when adsorbate-induced states are near the Fermi level, akin to observations in the metallic systems.Likewise, the sensitivity to σ occurs due to the different responses of slabs, with and without adsorbates, to changes in σ (Figure S6(a)).To minimize the smearing effect, extremely low σ values (≤ 0.005 eV) and the non-finite-temperature formalism were used (Figure S6(b)).However, it is important to note that such low σ values can lead to numerical instability.For the RPA methods prone to instability, such as RPA@PBE and RPA@HSE06, it is important to ensure numerical stability by verifying the consistency of DFT energy before and after the exact diagonalization step.Additionally, the RPA@PBE method requires denser k-point grids compared to those required by global hybrid-based approach, likely due to the influence of states near the Fermi level (Figure S6(b,c)).In RPA@PBEx approaches for non-metallic systems, the same or denser k-point grids, compared to the

Figure 1 :
Figure 1: Error of the calculated formation energies from the experimental formation energies.All formation energies are computed by linear regression to eliminate systematic errors.The grey shaded region represents chemical accuracy, ±0.0433 eV (±1 kcal/mol).PBE0-D3 plot is eclipsed by HSE06-D3 plot.

Figure 2 :
Figure 2: Error of the calculated reaction energies from the experimental reaction energies.Grey shade is a region of the chemical accuracy, ±0.043 eV (±1 kcal/mol).PBE0-D3 plot is eclipsed by HSE06-D3 plot.

Figure 3 :
Figure 3: Error of the calculated unitcell volume per atom from the experimental unitcell volume.

Figure 4 :
Figure 4: Error of the calculated adsorption energies in (top) metals, (center) metal oxides, and (bottom) MOFs from the adsorption energy of optRPA.The distribution of errors for each functional are plotted on the right as box and whisker plots.The top, middle, and bottom line of the box indicates 75% quartile (Q3), median, and 25% quartile (Q1) respectively.The whiskers extend from the edge of the box to the smallest and largest values within 1.5 times the interquartile range of Q1 and Q3, respectively.All three hybrid functionals and optRPA numbers are calculated by running single-point calculation on PBE-optimized geometry.Yellow stars in the plot indicate the reported experimental adsorption energies corrected to be directly comparable to DFT energies.For Pd, the adsorption energy from RPA@PBE is set to be the ground truth because optRPA calculations require significantly higher computational cost due to the use of the densest k-point grids.

Figure 5 :
Figure 5: Magnitude of the D3 correction as a function of the adsorption energy for all systems and functionals.Negative sign in the y-axis means that the adsorbate-adsorbent interaction is energetically favorable.

Figure S1 :
Figure S1: Error of the atomization energies from CCSD(T) values.

Figure S2 :
Figure S2: Error of the formation energies from CCSD(T) values.

Figure S3 :
Figure S3: Convergence test results of E c,RP A contribution to ∆ r E ads in terms of ECUT.*denotes the parameters used in the main text.

Figure S4 :
Figure S4: Convergence test of RPA@PBE calculations on Pd.(a) with E HF correction.(b) without E HF correction and after subtracting 1 2 T∆S elec .* denotes the parameters used in the main text.

Figure S5 :
Figure S5: Convergence test results of RPA@PBE0 calculations on Cu in terms of pseudopotentials (a) and k-points grid (b).* denotes the parameter used in the main text.

Figure S7 :
Figure S7: Error of calculated formation energy per each gas molecule without D3 correction.

Figure
Figure S8: Error of calculated reaction energy from experimental reaction energy without D3 correction.
Figure S8: Error of calculated reaction energy from experimental reaction energy without D3 correction.

Figure S11 :
Figure S11: Error of the calculated adsorption energies with (top) and without (bottom) D3 correction compared to the adsorption energy of optRPA in deprotonated MIL-125 and protonated MIL-125.

Figure S12 :
Figure S12: Relation between BEEF-vdW error bar and standard deviation of adsorption energies from different functionals.Rgression line is drawn, and confidence interval of 95% is shaded in red.Pearson correlation coefficient is 0.536.

Figure S13 :Figure S14 :
Figure S13: Relation between adsorption energy and D3 correction with regression lines drawn.Confidence interval of 95% is shaded in color.

Table 2 :
Max error in magnitude, mean signed error, and root mean squared error of the calculated unitcell volume per atom from the experimental unitcell volume.Numbers in paranthesis are the errors calculated without MIL-125 since thermal corrections are not available for MIL-125.OCUPUY is not included in any error calculations due to inconsistent geometries.Units are Å3 /atom.

Table 3 :
Error of the calculated adsorption energy compared to the adsorption energy of optRPA.Unit in eV.

Table 4 :
Standard deviation and error of the calculated adsorption energy compared to the adsorption energy of optRPA, for each material."std.dev" stands for standard deviation between functionals (except for RPA results).Unit in eV.

Table 5 :
Comparison of discrepancy in calculated adsorption energies between geometry-constrained OCUPUY (fixed PBE geometry) and unconstrained OCUPUY (original), relative to fully constrained OCUPUY (fixed PBE geometry and PBE magnetism).Unit in eV.

Table S1 :
Three sets of different kspacing and kinetic energy cutoff used for testing convergence of energy

Table S3 :
Convergence test result of bulk materials' energy in eV/atom unit.One GGA (PBE-D3), one meta-GGA (SCAN-D3) and one hybrid functional (PBE0-D3) are chosen for the test.The most accurate setup (A3) is set to be zero, so the numbers to be the energy difference from the A3 setup.The convergence criterion is 0.025 eV/atom.

Table S4 :
Convergence test result of gas molecules' energy in eV unit.One GGA (PBE-D3), one meta-GGA (SCAN-D3) and one hybrid functional (PBE0-D3) are chosen for the test.Only Γ point is sammpled, and the kinetic energy cutoff of 400 eV, 500 eV, 600 eV, 700 eV are tested.The most accurate setup (700 eV) is set to be zero, so the numbers to be the energy difference.The convergence criterion is 0.025 eV, but to be aligned with the convergence test of bulk materials, 600 eV is chosen for the calculations of gas formation energy and gas reaction energy.

Table S6 :
Convergence test results of RPA molecular atomization energy in eV unit.The most accurate setup (700 eV and 600 eV for E HF and E c,RP , respectively) is set to be zero, so the numbers to be the energy difference.

Table S7 :
Error of the calculated formation energies and reaction energies compared to the experimental formation and reaction energies without D3 correction.∆LS f E − ∆ f E exp ∆ r E − ∆ r E exp

Table S8 :
Error of calculated unitcell volume per atom from the corrected experimental volume with and without D3 correction.Unit in Å3 /atom.V DF T − V exp with D3 V DF T − V exp without D3

Table S12 :
Standard deviation and error of the calculated adsorption energy without D3 correction compared to the adsorption energy of optRPA, for each material."std.dev" stands for standard deviation between functionals (except for RPA results).

Table S13 :
Deviation of calculated adsorption energies from experimental values (in eV), with and without D3 correction.

Table S14 :
Comparison of adsorption energy and interaction energy (same gas and solid geometry from the RPBE optimized gas+solid system, and compute single-point calculation to get energies), and the amount of D3 correction added.

Table S16 :
Results of different spin state configurations tested and their energy in OCUPUY, using PBE-D3.Numbers in 'Initial magnetic moments' and 'Optimized magnetic moments' represent magnetic moment of four vanadium atoms in the unit cell.3 represents −3.

Table S17 :
Different spin states and energy after gas adsorbed in OCUPUY, using PBE-D3.Original state represents the result of Figure4.'3333' state represents the state where all four vanadium atoms in OCUPUY unitcell have up-spin states.∆E stands for the energy difference of the '3333' state compared to the original state (unit in eV).'mag' stands for the total magnetic moment.Since NH 2 and CN are already in the '3333' state (or 3333 ), the last row of them is empty.