DFT Benchmark Studies on Representative Species and Poisons of Methane Steam Reforming on Ni(111)

Ni catalysts used in Methane Steam Reforming (MSR) are highly susceptible to poisoning by carbon-based species, which poses a major impediment to the productivity of industrial operations. These species encompass graphitic carbon-like formations that are typically modelled as graphene. First principles-based approaches, such as Density Functional Theory (DFT) calculations, can provide valuable insight into the mechanism of graphene growth in the MSR reaction. It is, however, critical that a DFT model of this reaction can accurately describe the interactions of Ni(111) with the MSR intermediates as well as graphene. In this work, a systematic benchmark study has been carried out to identify a suitable DFT functional for the graphene-MSR system. The binding energies of graphene and important MSR species, as well as the reaction energies of methane dissociation and carbon oxidation, were computed on Ni(111) using GGA functionals, DFT-D and van der Waals density functionals (vdW-DF). It is well-established that the GGA functionals are inadequate for describing graphene-Ni(111) interactions. In the case of vdW-DF, the optPBE-vdW functional predicts the binding energies of graphene and several important MSR species with reasonable accuracy; however, it provides poor estimates of CO and O binding energies. Among the DFT-D functionals, PBE-D3 and PBE-dDsC have been found to exhibit acceptable accuracy for graphene and most MSR species (excluding adsorbed CO), and therefore, both functionals are promising for elucidating carbon-based catalytic poisoning in the MSR reaction. Overall, no single DFT functional could estimate the binding energies of be around 0.29 ML. The reaction energy of “carbon oxidation” (to form carbon on Ni(111) and O 2 (g)) is derived using the experimental carbon adsorption energy (with respect to C(g)) and the CO gas-phase formation energy (obtained from the atomisation energies dataset of CCSD(T) 56 ). The value is reported at 0 K (the ZPE/thermal contributions of “carbon adsorption” reaction is removed using the PBE functional prediction).


Introduction
Methane steam reforming (MSR) is an important industrial process for the production of carbon monoxide and hydrogen (commonly referred to as syngas). These species are crucial in various industrial applications, for instance, in the Fischer-Tropsch process, syngas acts as a feedstock for manufacturing a variety of liquid hydrocarbons. 1 Hydrogen is widely used in energies of H and H 2 O on the Ni(111) surface. PBE was able to estimate the binding energy of hydrogen with acceptable accuracy, however, it significantly under-predicts the water binding energy. On the other hand, DFT-D3 functionals (such as PBE-D3 and RPBE-D3) predict the binding energy of water with excellent accuracy but overestimate the hydrogen binding energy. Göltl et al. 29 benchmarked the performance of important GGA functionals, DFT-D functionals and vdW-DF for predicting the heat of molecular/dissociative adsorptions of CH 4, CH 3 I, CH 3 , I and H on Ni(111). The authors conclude that no individual DFT functional has high accuracy for all the species. For instance, the PBE-D3 functional was found to exhibit high predictive accuracy for species such as CH 4 , CH 3 and CH 3 I, whereas the BEEF-vdW functional gave the best estimate for the binding energies of H and I. Gautier et al. 30 computed the binding energies of carbon monoxide, oxygen, hydrogen and several hydrocarbons on Pt(111) using DFT functionals such as PBE, PBE-dDsC and few other vdW-DF. The authors compared the predictions of these functionals with micro-calorimetry data obtained from the literature. In addition, Wellendorff et al. 31 carried out a benchmark study for the DFT functionals using experimental data of ten different transition metals. There are other DFT benchmark studies that have tested the applicability of DFT-D functionals and vdW-DF for various systems. 32,33 Nevertheless, to the best of our knowledge, a thorough investigation on the performance of GGA functionals, DFT-D functionals and vdW-DF in predicting the binding energies of key intermediates of MSR and graphene on Ni(111) is not available in the literature.
The main objective of this work is to do a systematic comparative study on the predictive accuracy of relevant GGA functionals, DFT-D functionals and vdW-DF for graphene and a few representative MSR elementary reactions. To achieve this goal, the binding energies of graphene and MSR species (CO, C, CH 3 , H 2 O, H, O and OH), and the reaction energies of methane dissociation (considered to be the rate-determining step of MSR) and carbon oxidation were computed on the Ni(111) surface. The DFT functionals include PBE 34 , RPBE 35 and revPBE 36 (GGA functionals), PBE-D3, RPBE-D3, revPBE-D3 16 , PBE-dDsC 18,19 and PBE-TS 17 (DFT-D functionals), optB86b-vdW, optB88-vdW, optPBE-vdW 22 and BEEF-vdW 24 (vdW-DF). The predictions of the aforementioned DFT functionals have been benchmarked against experimental and computational data (from higher-level theory) available in the literature.
The rest of this manuscript is organised as follows. The "Methods" section includes the computational details of the DFT calculations and the relevant formulae used to estimate the DFT binding energy, the zero-point energy (ZPE) correction, the thermal energy contributions and the root mean-square deviation (RMSD). Subsequently, the "Results and Discussion" section presents the comparison of DFT predictions of graphene and MSR species to relevant literature data and puts the key findings in the context of practical catalysis. Finally, in the "Conclusions" section, a broad perspective on the performance of GGA functionals, DFT-D functionals and vdW-DF is provided along with some needs and potential opportunities for DFT method development are also discussed.

Methods
Spin-polarised DFT calculations have been performed using the Vienna ab initio Simulation Package (VASP). The projector augmented wave (PAW) method was used to model the interactions between core and valence electrons. A plane-wave basis set was employed, and the kinetic energy cut-off was set to 400 eV (refer to Figures S1(b) and S2(b) of the supplementary information (SI) for the plane-wave cut-off convergence plots). The Ni bulk calculations were carried out using a 191919 k-point grid. In the Ni bulk system, the electron smearing is performed using the tetrahedron method with Blöchl corrections and the smearing width was set to 0.05 eV. The Ni(111) surface was modelled with a six-layer p(44) slab, where the Ni atoms of the three bottom-most layers were fixed at the optimised lattice constant and others were allowed to fully relax. The periodic images along the z-direction were separated by a vacuum of 12 Å. The first Brillouin zone of Ni(111) was sampled with a 551 Monkhorst-Pack k-point grid. The k-point mesh was chosen by performing convergence tests on carbon-Ni(111) and graphene-Ni(111), respectively (refer to Figures S1(a) and S2(a) of the SI). The electronic self-consistency tolerance was set to 10 -7 eV, and the geometry optimisation terminated when the Hellmann-Feynman forces acting on the atoms that are allowed to relax reached a value less than 0.01 eV/Å. The smearing of electrons was carried out using the Methfessel-Paxton method and the smearing width was set to 0.1 eV. In the vibrational calculations, the Hessian matrix was computed using the central finite-difference method. The step size of the displacement was set to 0.02 Å. (1) = + ( ) - When an adsorbate binds to the Ni(111) surface, the translational and rotational motions are reduced to vibrational modes. The thermal energy contribution of the bound state has been computed by using the harmonic approximation. 37,38 According to Réocreux et al. 39 , the thermal corrections of chemisorbed species obtained under the harmonic approximation have acceptable accuracy (even for low-frequency modes). The bound state vibrational calculation is performed by fixing the positions of the Ni surface atoms (vibrational modes of the Ni surface atoms, commonly referred to as phonons, are assumed to have an insignificant contribution to energies of adsorption or reaction). Equations 5 and 6 provide the thermal correction for species in the gas phase and the bound state, respectively. In equation 7, the ZPE and thermal corrections are added to the DFT predicted binding energy.
The term from equation (7) is used to make comparisons with experimental data and ℎ the root mean square deviation (RMSD) defined by equation (8)  represents the experimental binding energy of the i th species.

Results and Discussion
As discussed in the "Introduction section", the whisker carbon poisoning of MSR primarily involves the formation of graphene islands on the support side of the

DFT benchmarks of graphene-Ni(111)
For all of our calculations of pure graphene and Ni(111), we have used the DFT-predicted lattice constants which are provided in the SI for each functional (refer to Table S1). For the calculations of graphene bound to Ni(111), the Ni lattice constant was adopted. It is known that the LDA, GGA and hybrid functionals fail to capture the van der Waals interactions between graphene and Ni(111) 15 , and thus, in this benchmark study, the performance of DFT-D and vdW-DF has been tested.
The DFT predictions must be compared to reliable experimental or theoretical values; however, there is limited experimental data on the binding energy of graphene on Ni(111).
Shelton et al. 40 have studied the segregation behaviour of carbon on Ni(111) and concluded that the binding energy of graphene on Ni(111) is around 50 meV greater than the graphite exfoliation energy. Recent experimental works have reported the exfoliation energy of graphite to be in the range of 24.87 meV/C atom -66.33 meV/C atom. 26,41,42 Based on the findings of Shelton et al. 40 , the graphene binding energy on Ni(111) can be deduced to be in the range of 74.62 -116.08 meV/C atom. 26 As the graphene binding energy estimate is derived using the findings of only a single experimental study (Shelton et al. 40 ), a certain level of caution must be exercised while comparing this value to the DFT predictions. Thus, we also compare our calculated values to interactions. 14,43 The RPA prediction of graphene binding energy will therefore be used in assessing the performance of DFT functionals. The RPA predicts two minima for graphene on Ni(111), located at graphene-Ni distances of 2.17 Å and 3.25 Å, respectively. 14,44,45 The RPA prediction at the former minimum (referred to as "the first minimum" henceforth) is of importance as it is in close agreement with the experimental graphene-Ni distance of 2.11 ± 0.07 Å. 46 Figure 1  configurations). 47,48 In Table 1, the DFT binding energy predictions of this work are compared to the values reported by other theoretical studies available in the literature. The graphene-Ni distance predictions of the DFT functionals are also recorded in Table 1. 1). Furthermore, Mittendorfer et al. 14 have shown that the optB88-vdW functional predicts the binding energy of graphene with high accuracy. As illustrated in Figure 1, the optB86b-vdW functional substantially overpredicts the binding energy of graphene with respect to the RPA value, though the predicted value is within the experimental range. Shepard and Smeu 44 and Li et al. 25 have reported similar binding energy values for graphene using the optB86b-vdW functional. On the contrary, the calculations performed by Janthon et al. 26 show that there is close agreement between the optB86b-vdW prediction and the RPA value (refer to Table 1 for the optB86b-vdW predictions obtained from the literature).

Table 1. Comparative study of DFT predictions of the first minimum of graphene-Ni(111)
Ref. 14 uses a five-layer Ni (111) Table   S16 and Figure S7).
As shown in Table 1, the optPBE-vdW functional estimates a binding energy value of 41.79 meV/C atom for graphene-Ni(111), which is in reasonable agreement with the RPA value.
Janthon et al. 26 report that the optPBE-vdW functional predicts a much weaker adsorption energy for graphene-Ni(111) (9.33 meV/C atom setup compared to our study; this could be the reason for such a substantial difference in the reported optPBE-vdW predictions (key information on the computational setup of the literature studies is provided in the footnote of Table 1). The BEEF-vdW functional was found to predict a weak repulsive interaction between graphene and Ni(111), as also pointed out by Wellendorff et al. 24 The computational error estimate of BEFF-vdW is large (obtained by employing an ensemble of 2000 functionals 24 ). Certain caution needs to be exercised while using the BEEF-vdW error estimate which results in a graphene binding energy range from -76.54 meV/C (negative values are qualitatively incorrect with respect to the RPA prediction) to 74.36 meV/C (substantially strong binding affinity between graphene and Ni(111); which closely matches the RPA value). This wide variation also underscores the importance of using the appropriate exchange-correlation approximation to obtain accurate results.
The PBE-TS functional significantly overestimates the graphene binding propensity on the Ni(111) crystal (refer to Figure 1). As shown in Table 1    Taken together, the results of our benchmark studies show that the graphene-Ni (111) interactions are best represented by the optB88-vdW functional. The PBE-D3, PBE-dDsC, and optPBE-vdW functionals also perform reasonably well.

DFT benchmarks of MSR species
In a second set of benchmark calculations, the DFT functionals were used to estimate the reaction energies of methane dissociative adsorption, carbon oxidation, and the binding energies of representative species of the MSR reaction, in particular: CO, C, CH 3 , H 2 O, H, O and OH on Ni(111). The predictions of these functionals have been systematically compared against experimental data, the sources of which are compiled in Table 2. The binding sites of the MSR species are provided in the SI (refer to Figure S4). For each MSR species under study, Table 2 reports the pertinent adsorption reaction, the details of the experiment that we compare against, the simulation setup of the corresponding DFT calculation (this information is provided as a footnote in Table 2), and finally, the value of the experimental binding energy.
In Table 2, ML is defined to be the ratio of number of adsorbate atoms/molecules to the number of Ni atoms on the Ni(111) surface. An important caveat to be noted is that the experimental surface coverage of carbon is not clearly known (refer to the last row of Table 2 for more details), therefore, certain caution is exercised while comparing DFT predictions to the experimental binding energy of carbon.   b) The experimental internal energy of adsorption is reported (the heat of adsorption obtained from the experiment is reduced by k B T). 57 This value is compared to the DFT prediction . ℎ c) The Arrhenius activation barrier obtained from flash desorption studies of hydrogen was compared to the ZPE corrected DFT energy prediction (no thermal corrections were included, only equations 2-4 of the "Methods section" were used in this case) d) Since the experimental surface coverage is 0.5 ML, a p(22) Ni(111) supercell (constituting two OH atoms adsorbed on the three-fold hollow sites) has been used and the Brillouin zone integration was done with a 10101 k-point grid. e) The kinetic isotope effects are assumed to be negligible. f) By convention, a positive value indicates exothermic reaction, in which case the dissociative adsorption is favourable. A negative value indicates endothermic dissociative adsorption (not favoured).

System
Experimental method (eV) CO(g) + *   CO* The differential heat of adsorption of CO on Ni(111) is recorded at the limit of zero-coverage using calorimetry. The temperature is maintained at 300K.
The associative desorption of hydrogen is studied using flash desorption spectroscopy. The desorption of H/Ni(111) has been observed to follow second-order kinetics at low coverages. The Arrhenius activation barrier of associative desorption of hydrogen has been estimated.
The reaction energy of "carbon oxidation" (to form carbon on Ni(111) and O 2 (g)) is derived using the experimental carbon adsorption energy (with respect to C(g)) and the CO gas-phase formation energy (obtained from the atomisation energies dataset of CCSD(T) 56 ). The value is reported at 0 K (the ZPE/thermal contributions of "carbon adsorption" reaction is removed using the PBE functional prediction). In Table 3, the experimental values and the DFT predictions of MSR species are systematically recorded. The colour code in Table 3 gives an indication on the extent of deviation of the DFT prediction from the corresponding experimental value (the convention of the colour code is stated as a footnote in Table 3). As shown in Figure 4, the deviations of the DFT predictions from the experimental values have been presented using radar plots.
Furthermore, the RMSD values have been computed to make a quantitative assessment on the overall performance of the DFT functionals (refer to Figure 5).
The oxygen binding energy predictions (reported in  Table S3). Similarly, the O 2 (g) DFT energy has been used to calculate the reaction energy of carbon oxidation (C* + (½)O 2 (g)   CO(g)).
The PBE functional predicts the reaction energy of methane dissociation with reasonable accuracy (the deviation is -0.23 eV). It performs appreciably well for MSR species such as hydrogen, oxygen and hydroxyl. The PBE functional overpredicts the binding energy of carbon monoxide by more than 0.5 eV (as depicted in Figure 4(a)). In the literature, it is well established that the PBE exchange-correlation approximation fails to accurately predict the binding energy of carbon monoxide. 31 The PBE functional underpredicts the binding energy of H 2 O, as it does not include the dispersive interactions between H 2 O and Ni(111). Zhu et al. 28 also report the same behaviour for the PBE functional. As shown in Table 3 and Figure 4(a), the PBE functional accurately estimates the carbon binding energy and the reaction energy of carbon oxidation (however, we note again that the experimental value of carbon might be subject to finite coverage effects).
The RPBE and revPBE functionals give a better estimate of CO binding energy than the other DFT functionals tested in this study (as shown in Figure 4(a) and  Table 3). As shown in Figure 5, the PBE functional has better overall performance than other GGA functionals. Nonetheless, the GGA functionals do not account for the dispersive interactions of graphene-Ni(111) (as stated in the "Introduction" section), and therefore, they are not suitable to investigate the carbon poisoning mechanism of Ni in the MSR reaction.

Table 3. Comparisons of DFT predictions with experimental values of MSR species
Deviation = |DFT predicted value -Experimental value|. The colour code has the following convention: green:  It can be inferred from Table 3 and Figure 4(b) that the performance of functionals within the vdW-DF class is unsatisfactory. Although optB86b-vdW and optB88-vdW functionals predict the reaction energies of methane dissociation and carbon oxidation, and the binding energies of CH 3 , H 2 O, H and C with acceptable accuracy, both functionals perform poorly for species such as CO and O (as depicted in Figure 4(b)). In agreement with the findings of our study, Hensley et al. 58 also report that the optB86b-vdW and optB88-vdW functionals do not predict the binding energies of CO with acceptable accuracy. The RMSD values of optB86b-vdW and optB88-vdW functionals are 1.11 eV and 1.00 eV, respectively, which is much higher than most of the other dispersion-inclusive DFT functionals (refer to Figure 5). Thus, neither of these functionals appear promising for studying graphene formation in the MSR reaction.  Table 3). The BEEF-vdW functional has the highest RMSD score among all the dispersion-inclusive DFT functionals benchmarked in this study (refer to Figure 5); which makes it unsuitable for the graphene-MSR system. The optPBE-vdW functional predicts the binding energy of OH with high accuracy. For species such as H 2 O, H and C, the optPBE-vdW functional produces deviations within an acceptable range (0.14-0.16 eV). It estimates the reaction energy of methane dissociation with reasonable accuracy (the deviation is 0.19 eV). The carbon oxidation reaction energy and the CH 3 binding energy are under-predicted by 0.11 eV and 0.27 eV, respectively. In the case of species such as CO and O, the predictions of optPBE-vdW are in poor agreement with experimental data (refer to Table 3 and Figure   4(b)). Nevertheless, the optPBE-vdW functional has a much better overall performance. The RMSD score of optPBE-vdW functional is 0.85 eV, which is substantially lower than the other vdW-DF (as shown in Figure 5). Furthermore, as discussed in section 3.1, the optPBE-vdW functional predicts the first minimum of graphene with reasonable accuracy. Therefore, among the vdW-DF, the optPBE-vdW functional can be considered to study the chemistry of carbon poisoning on Ni(111) (with appropriate corrections to CO and O binding energy values).
In contrast to vdW-DF, the DFT-D functionals perform reasonably well for this system. The The PBE-D3 functional predicts the binding energies of OH and C with appreciable accuracy.
It provides a good description of the H 2 O-Ni(111) interactions, whereas, the hydrogen binding energy is substantially overestimated; Zhu et al. 28 have also made similar observations. According to Table 3, PBE-D3 predicts the CH 3 binding energy to be 2.25 eV, which agrees closely with the experimental result. Furthermore, it has been found to predict with reasonable accuracy the reaction energies of methane dissociation and carbon oxidation, as well as the oxygen binding energy (the deviations are less than 0.25 eV; refer to Figure 4(c)). In a recent study, Göltl et al. 29 Figure   4(c) and Table 3).  Figure 5). They also account for the graphene-Ni (111)  promising for more detailed studies to understand the carbon poisoning chemistry of Ni (111) in the MSR reaction (suitable corrections are required for the CO binding energy predictions of these functionals).

Conclusions
Detailed DFT studies of the MSR reaction on Ni can provide valuable mechanistic insight on catalyst poisoning and deactivation phenomena, which have a severe impact on industrial operations. At the molecular scale, the poisoning mechanism primarily involves the formation of graphene islands; thus, the accurate description of graphene-Ni (111)