Hydrogen adsorption on transition metal carbides: a DFT study.

Transition metal carbides are a class of materials widely known for both their interesting physical properties and catalytic activity. In this work, we have used plane-wave DFT methods to study the interaction with increasing amounts of molecular hydrogen on the low-index surfaces of four major carbides - TiC, VC, ZrC and NbC. Adsorption is found to be generally exothermic and occurs predominantly on the surface carbon atoms. We identify trends over the carbides and their surfaces for the energetics of the adsorption as a function of their electronic and geometrical characteristics. An ab initio thermodynamics formalism is used to study the properties of the slabs as the hydrogen coverage is increased.


Introduction
One of the key scientific and technical challenges of the 21st century is the need to meet the global demand for energy in a clean, secure and sustainable way, 1,2 by overcoming remaining issues of electricity storing and pollution. 3 The hydrogen evolution reaction (HER) has been intensively studied as a route to produce H 2 from H 2 O, with the primary goal of using molecular hydrogen in carbon-free fuel cells and as a reductant to recycle CO 2 . The HER is most effectively catalysed by Pt electrodes, which allow the reaction to be driven using a low overpotential. 4,5 However, the high cost and low supply of Pt make it an impractical catalyst for global scale production of H 2 .
Low cost materials have been investigated in recent years with the aim of obtaining similar performances to Pt. 6,7 One such avenue of research is the class of materials known as transition metal carbides (TMCs), which couple low cost and high durability with interesting catalytic properties for many reactions, [8][9][10][11] as originally highlighted by the seminal work by Levy and Boudart in 1973. 12 Transition metals can form carbides in a variety of structures and stoichiometries, depending on their position in the periodic table and the extent of defects in the structure. 13,14 Monocarbides and non-stoichiometric carbides of early transition metals are found to have stable rock-salt bulk structures, [15][16][17][18][19][20] while metals at the centre of the transition series have more complex stoichiometries and structures. 11,21,22 These compounds combine some of the physical properties of transition metals, ionic and covalent solids, 23,24 which can be further modified by tuning the amount of carbon in the lattice and the bonding character of the metal-carbon bond. 25,26 The catalytic behaviour of TMCs has been investigated for a wide range of reactions, such as hydrogenation, CO 2 reduction, water-gas shift, hydro-treatment and a number of processes involving hydrocarbons, including isomerization and methane reforming. [27][28][29][30][31][32][33] For the hydrogen evolution reaction, these materials have been used with promising results as both supports for platinum electrocatalysts 33 and as catalysts themselves, either as bulk materials or nanoparticles. 11,21,34 However, TMC catalysts have still not seen widespread adoption in industrial processes. 35 In this work, the interaction between molecular hydrogen and four TMCs from groups IV and V is studied using planewave Density Functional Theory (DFT) methods. TiC, ZrC, VC and NbC, all with 1 : 1 stoichiometry and the rock-salt structure, are modelled and cut into four different low index surfaces, i.e. the (001) and (011) surfaces and the carbon-and metalterminations of the (111) plane. We have identified trends in the energetic and electronic effects of the H 2 adsorption from which we can make predictions concerning the H 2 evolution under specific conditions on the different surfaces. The metal carbides were chosen after a comprehensive review of bulk and surface properties showed that they had promising properties as potential catalysts for CO 2 hydrogenation reactions. Previous experimental and computational results 27,36 suggest that TiC and its immediate neighbours are good candidates for H 2 evolution and CO 2 reduction and could be superior to other transition metal carbides with 1 : 1 stoichiometries. Whilst these prior studies have clearly provided important insight, our work offers the first comprehensive theoretical study into hydrogen adsorption on all low-index surfaces of TiC, ZrC, VC and NbC.

Computational details
All calculations have been performed within the framework of the density functional theory (DFT) as implemented in the VASP (Vienna Ab initio Software Package) code. [37][38][39] The electronic energy was obtained using PAW (projected augmented wave) potentials 40 and plane-wave basis sets, respectively, for core and valence electrons, using the Perdew-Burke-Ernzerhof (PBE) functional. 41 Corrections were added to account for long-range interactions by the semi-empirical Grimme D3 dispersion method 42 and for non-spherical contributions to the PAW potentials as incorporated in the code. All energies were converged within a cutoff of 520 eV and an electronic self-consistent field (SCF) threshold of 10 À5 eV. Convergence was determined using the tetrahedron method, implementing Blochl corrected smearing 43 and in all cases spin polarization was enabled, with the only exception of the isolated H 2 closed shell molecule.
The surfaces of the TMCs in the study were simulated by 2 Â 2 Â 3 supercell slab models generated via the METADISE code, cutting the bulk along the (001), (011) and (111) planes. Each slab has 6 atomic layers and 16 atoms per layer. The (001) and (011) surfaces are created such that they preserve the bulk stoichiometry, resulting in an equal number of carbon and metal atoms being exposed to the vacuum. Conversely, the (111) plane is parallel to carbon and metal layers, resulting in two possible surface terminations that respectively expose carbon and metal atoms in the surface. Such a protocol has been applied previously to TMCs and has been shown to describe accurately the electronic and structural properties of the (111) surfaces, as will be further discussed below. 30,36 To avoid interaction along the axis perpendicular to the surface, 12 Å of vacuum were added in that direction. The 5 Â 5 Â 1 K-points reciprocal lattice matrix was generated using the Monkhorst-Pack method. The energy of the system was minimized by keeping the cell parameters fixed at their bulk-optimised value and allowing relaxation of all atoms, apart from those belonging to the 2 bottom layers of the TMC, which were kept fixed. The minimum energy structures were found using a built-in DIIS algorithm with a convergence force threshold of 10 À2 eV Å À1 .
The interactions between the different surfaces and hydrogen were compared using the adsorption energy (DE ads ) as calculated by: where E slab+nH is the energy of the hydrogenated slab, E slab and E H 2 are the energies of the pristine surface and isolated H 2 , respectively, and n is the number of hydrogen atoms involved in the adsorption. By normalizing, we consider as DE ads the adsorption energy per atom instead of the total energy change in the system. The energy of hydrogen, which is also used to calculate the chemical potential of the molecule, is calculated in a broken symmetry unit cell of 9 Â 10 Â 11 Å, using the same parameters with which the surfaces have been investigated. Atomic charges were calculated on specific configurations through a Bader charge analysis, 44 as implemented by Henkelman and co-workers. 45,46 The reported work functions (F) were calculated by subtracting the Fermi energy from the vacuum energy, extracted from the local potential by averaging the energy levels above the surface. The resulting F were further used to calculate the theoretical half-cell electrode potential for the slabs, for each coverage state, as outlined in the seminal paper by Trasatti: 47 where F slab is the work function of the slab, e is the electron charge and 4.5 V is an estimate of the absolute potential of the NHE as calculated from its work function, coherently with that found in the literature. 48,49 The effect of hydrogen coverage has been further investigated by means of the ab initio thermodynamics formalism 50,51 with the objective of elucidating the level of hydrogenation as a function of the chemical potential of H 2 in the gas phase. The surface free energies (s) are therefore defined as follows: where g is the surface energy of the relaxed, pristine slab, E y is the energy of the slab with coverage state y, E clean is the energy for the pristine slab, n is the number of H atoms adsorbed on the surface and m H is their chemical potential in the gas phase.
Values for the chemical potential as a function of temperature and pressure were obtained using standard statistical thermodynamics, via the formula: where E el and ZPE are the electronic and zero point energy of the hydrogen molecule as derived from DFT, the DG term is the difference in Gibbs free energy between 0 K and T at constant pressure p 0 = 1 bar as calculated from the partition function of the molecule and the last term is the pressure dependency of the Gibbs free energy.

Results and discussion
Our study of the H 2 adsorption on these TMCs -TiC, VC, ZrC and NbC -first explored the stability of molecular hydrogen above a test surface, namely the (001) surface of TiC. Initially, H 2 was placed 6.5 Å over the (001) surface of TiC and fully optimized, resulting in a weak long-range interaction energy of 0.05 eV. Secondly, the distance between the surface and the molecule was systematically decreased to find potential curves. Spontaneous dissociative chemisorption was observed to occur when the molecule was placed less than 2.0 Å above the topmost surface atomic plane. The dissociation happens on top of one of the surface Ti atoms, producing an intermediate in which one of hydrogens is coordinated by two Ti atoms and the other absorbs on top of a C atom. However, the most favourable adsorption is found when both H atoms absorb on top of a surface C, a feature that, as shown below, is shared by all tested systems. The same protocol was subsequently repeated on

Hydrogen adsorption
Having established the fundamental exothermicity of adsorption, adsorption sites and their respective energies were explored on all carbides. The results show clear geometric and energetical similarities between each surface plane amongst the four studied carbides. Two different adsorption geometries for the (001) surface are shown in Fig. 1; in the first the adsorbed hydrogens are located on top of a C atom (C on top site) and in the second between two metal atoms (M 2-hollow site). On TiC, in the C on top site, each H atom absorbs 1.28 Å above the surface, while on the M 2-hollow site, the H absorbs directly above the surface Ti. In the first case, the formation of a C-H bond is combined with a change in M-C bond distance, which increases slightly, as the dangling bonds of C are linked with the electron on hydrogen. However, this elongation is very small. Energetically, the C on top site is favoured, with DE ads = À0.62 eV, while the M 2-hollow site shows an endothermic adsorption, with DE ads = +0.38 eV. A Bader charge analysis and a density of states (DOS) calculation have been performed on both the pristine and the hydrogenated surfaces. The two analyses show how the adsorption has little effect on the electronic structure of the surface, and the former also highlights that for both adsorption sites the two adsorbed hydrogens retain one electron each. For further details on the two analyses, see ESI. † Similar results in terms of geometry, surface reconstruction and charge localization were found on all other carbides, as reported in Table 1. Energetically, TiC, VC and ZrC exhibit favourable hydrogen adsorption on their (001) surfaces, while NbC showing endothermic adsorption energies on both the C on top and M 2-hollow adsorption sites.
On the (011) surfaces, similar M 2-hollow and C on top adsorption sites are found. On such surfaces both sites show higher DE ads by 0.5-1 eV compared to (001), and the adsorption on the C on top site is again the most exergonic and thus most favourable on these surfaces. Due to the conformation of the surface, the C-H distance for the C on top site surface is shorter than that of the same site on (001) and an even smaller M-C bond elongation is found, while on the M 2-hollow site adsorbed hydrogen is at a greater distance from the surface. The distortion of the surface also causes the hydrogen atoms on the C on top site not to be vertically aligned with their respective carbon. The Bader charge analysis shows an equal localization of the charge between the two atoms, which as in the previous case retain their respective electrons. Table 2 shows the adsorption energies, Bader charges on hydrogens and C-H bond distances on the (011) surfaces of all four carbides.
The two surface terminations along the (111) plane are treated separately because of the variation in the stoichiometries of their topmost atomic layer and their resulting different surface energies. Both the metal-and carbon-terminated slab models show a dipole perpendicular to the surface, arising from the presence of alternate layers of carbon and metal atoms. 54 However, previous studies have shown similar trends in the electronic properties of the surface on both reconstructed and unreconstructed surfaces. 36 Moreover, experimental reports only show unreconstructed metal-terminated surfaces for the (111) plane, when this surface can be shown to be clean and free of impurities under ultravacuum conditions. [55][56][57][58] For these reasons, only the unreconstructed metal and carbon terminations of the (111) plane will be considered in this study.
On the metal-terminated (111) plane -(111)-M from now on -the surface atoms are arranged in a hexagonal fashion, as shown in Fig. 2. A threefold site is available for the hydrogen,    Table 3.
A different picture emerges from the calculations on the (111)-C surface. This surface termination shows instability in the presence of hydrogen, leading to a dramatic reconstruction on TiC and VC, whereby a combination of carbon atom sintering and hydro-carbon formation distorts the pristine surface geometries. However, for the two 5th period TMCs in this study, it was possible to find a minimum in the potential energy surface associated with the adsorption. In the latter cases the adsorbed hydrogen absorbs on a C on top site, similarly to that found on the (001) and (011) surfaces. However, a more exothermic reaction and a more extensive surface reconstruction of the model slab are involved. Table 4 summarizes these results, which are consistent with the previously reported energetic preference for metal termination on the (111) surface of these two carbides, as deduced from surface energies. 36 The binding energies obtained for hydrogen adsorption can be compared with previous studies on the same surfaces, 36 highlighting correlations between the surface properties and their reactivity. The reported surface energies and d-band positions for these materials show trends along the rows and down the groups that correlate well to their reactivity. The trends over the surfaces for a single carbide and over the same face of each carbide show an inverse correlation with the DE ads which, despite some discrepancies, can be considered a good descriptor of the activity of the carbide surface considered. Fig. 3 shows how, for the (001) surfaces of all carbides as well as all four surfaces of ZrC, a less negative d-band centre corresponds to a more exothermic hydrogen adsorption.

Higher hydrogen loadings
In this section, our approach is to load hydrogen atoms up to monolayer (ML) coverage. In the previous section, 16 adsorption sites had been identified on each model slab: on the slab representing (001) and (011) surfaces these are equally divided between C on top and M 2-hollow , while on those representing the two terminations of the (111) surfaces either 16 C on top or 16 M 3-hollow sites are present. Therefore, loadings of 4, 8 and 16 hydrogen atoms over the 16 surface adsorption sites are investigated, and the corresponding coverage and surface properties are then evaluated for each system.
On the (001) surface termination, adsorption of 4 and 8 H produces a decrease in adsorption energy per atom on each of    the four carbides in study. The adsorbed hydrogens all absorb on the 8 C on top sites. Subsequent adsorption of H 2 on the surface is expected to occur on the M 2-hollow site, which has been shown previously to be endothermic. However, steric effects prevent these sites from being available. Thus we find that on the carbides of the 5th period, ZrC and NbC, hydrogen adsorbs on the more endothermic M on top sites, resulting in a positive total adsorption energy (see Table 5 and Fig. 4). On the carbides of the 4th period, TiC and VC, it was not possible to find a stable minimum for the adsorption on such sites for the (001) surface. For these reasons we consider that the chemisorption of H on C on top defines the monolayer on all (001) surfaces in our study, even if these are only half of the stable minima found on the potential energy surface for chemisorbed hydrogen. Within our model slabs, this translates in having y = 1 when 8 H atoms are adsorbed.
For all four carbides, the increase in coverage leads to a decrease in the adsorption energy per atom, which is consistent with the results found for other carbides by previous work. 59 This effect is not unexpected, and insight comes from the analysis of the work function (F), which steadily decreases with coverage as electrons are inserted in the conduction band (Fig. 4). As highlighted in previous literature, 49,60 F can be correlated with the open circuit half-cell potential. We therefore used eqn (2) to convert the F values into potential values against the NHE. Such a conversion highlights how the potential decreases slightly with coverage. In particular, the voltage difference between pristine surface and the first adsorbed molecule is only 0.13 V, suggesting by reversing the relationship, how a small change in potential might influence the adsorption and desorption of H 2 on these surfaces.
The trends in adsorption energy, work function and theoretical V half cell with y are found to be similar on all the investigated carbides and surfaces. As shown in Table 5, adsorption up to the monolayer is exothermic on all investigated carbides which showed exothermic adsorption for a single molecule.
To elucidate the relationship between the different coverage conditions we computed the surface free energies of the hydrogenated surfaces. Fig. 5b shows the surface free energies for the four stable coverages of the TiC(001) surface versus the chemical potential of the hydrogen, whose dependency on temperature and pressure is shown in Fig. 5a. The graph highlights how lower coverages are expected to be more favourable for low chemical potentials. The vertical black line on the graph refers to the conditions of standard temperature and pressure, highlighting how these correspond to a fully hydrogenated surface, a property shared with every other investigated surface apart from the (001) surface of NbC, shown in Fig. 5c: on this surface, the increase in electronic energy as a result of the hydrogen adsorption results in a clean surface in standard conditions. However, both partial and full coverage are observed at higher chemical potentials, occurring at high pressures or low temperatures.
It is possible to correlate the catalytic activity of a material towards the hydrogen evolution reaction (HER) with the free energy of adsorption of hydrogen on its surface, which has to be close to zero to achieve the best catalytic activity. 61 This is because the reaction is composed of two separate steps, of which the adsorbed hydrogen is the necessary intermediate: at equilibrium the initial and final step will have the same free energy. Therefore, we can predict the conditions at which the reaction will be most efficiently catalysed by identifying the intersections between the lines representing two subsequent coverage states. On TiC(001) this would imply having the best catalytic activity at temperatures higher than 300 K, while on NbC(001) low temperatures or high pressures are needed to catalyse the reaction efficiently. As highlighted above, the graphs for VC(001) and ZrC(001) are similar to that of TiC(001), but due to the less exothermic DE ads values, the equilibrium states between different coverage states for VC(001) are closer to 298 K and 1 atm.
Similar trends are found on the (011) surfaces as a function of y, showing a decrease in DE ads , F and V half cell , as shown in Table 6. However, owing to the larger interatomic spacing on the (011) surfaces, the 8 M 2-hollow sites are still available after all 8 C on top sites are filled. As seen for the loading of a single molecule, the adsorption is exothermic on these sites, although less than the C on top sites.
The condition of full coverage on the studied (011) surfaces can therefore entail the filling of all metal and carbon sites, corresponding to 16 chemisorbed H atoms on 16 available adsorption sites in our model slabs. The calculated surface free energies show that such coverage is the most favourable for all but the most negative chemical potentials, corresponding to very high temperatures and low pressures (as shown in Fig. 5d for VC). We can therefore predict that these will be the conditions at which the (011) surfaces of the four carbides in our study will show HER activity, since less extreme conditions show strong bias towards the highest possible coverage. This effect is due to the larger DE ads found on the (011) surfaces, which make the adsorption favourable in all but the most extreme cases. The (111) surfaces have 16 equal adsorption sites on each slab, either C on top or M 3-hollow depending on surface termination.
On the metal termination, the adsorption is repeated without changes in geometry as hydrogen loading increases, which is followed by the expected decrease in DE ads , as shown in Table 7. On the carbon termination, higher loadings of hydrogen markedly stabilize the surface, revealing stable configurations for TiC and VC for 16/16 and 4/16 occupied adsorption sites, respectively, while keeping the reaction very exothermic for all coverages and carbides, as highlighted in Table 7.
Calculating the free energies for these surfaces produces similar results for all the investigated metal carbides. Owing to the high DE ads , hydrogens fill the adsorption sites (coverage = 1 ML) on both terminations. This result means that, in the phase space we have explored, these carbides are not likely to favour hydrogen desorption and will therefore exhibit poor catalytic activity toward the HER reaction on their (111) surfaces. This is exemplified by Fig. 5e, which shows the corresponding plot for ZrC(111)-M.
The calculations on the hydrogenated surfaces in this study have shown a strong dependence of the energetics of the reaction on both the surface and the material. Overall, it was possible to identify several trends from the data analysed. The C on-top is the preferred adsorption site, where the process is always exothermic except for NbC(001). Metal sites are less favourable on all surfaces, but show a stronger interaction as the number of coordinated metals increases, so that the M 3-hollow sites allow a strong adsorption of the hydrogen on (111) surfaces. On both carbon and metal sites the energetic trends and adsorption properties on the four surfaces studied are consistent across the four carbides, showing very similar trends. For all surfaces and coverages, reaction on VC and NbC corresponds to less exothermic adsorption energies with respect to the same reactions over TiC and ZrC, showing a clear effect of the group in the periodic table.
The effect of the row is less marked, but in some cases the reduced lattice parameter for the carbides of the 4th row induces steric effects that hinder the adsorption. These trends could be in some cases related to the electronic characteristics of the surfaces. The low F of TiC and ZrC is in relation to their higher activities, as are the d-band centres for all surfaces, which have been proven to show a useful correlation with the activity of the material over the hydrogen adsorption and could therefore be used as descriptors. Similarly, the effect of coverage is clearly correlated with a drop in both DE ads and F: as more hydrogens are added on a surface, more electrons are added to the conduction band of the solid, thereby obstructing further hydrogenation of the slab. However, the highest coverage is consistently the most exothermic across all surfaces, apart from the already mentioned NbC(001), as the decrease in exothermicity with coverage is small compared to the total value of adsorption energy, especially for (011) and (111) surfaces. For this reason, such surfaces are predicted to show less catalytic activity towards the hydrogen evolution reaction than (001) surfaces which show less exothermic DE ads values.

Summary
We have presented a comprehensive DFT study into hydrogen adsorption on the four low index surfaces of TiC, VC, ZrC and NbC. The adsorption is shown to occur primarily on carbon sites and on the most highly coordinated metal sites on the surface. We have demonstrated that all carbides show the same trends across their four surfaces. For each carbide the most exothermic adsorption is found on the (111) surfaces; however, the carbon termination of this surface proved to be very unstable on TiC and VC. The only non-exothermic adsorption is found on the NbC(001) surface. On increasing the amount of hydrogen loading on each of the selected surfaces we found a corresponding decrease in adsorption energy per atom, linked with a similar decrease in the work function of the surface. Although the full coverage is always connected to a lower total energy, the actual coverage of the surface has been found to be generally depended on temperature and pressure of the environment. Using ab initio thermodynamics we were able to predict that the (001) surfaces of the four carbides in the study show the most promise as catalysts for the HER, as the change in coverage states with hydrogen chemical potential is linked to thermodynamic equilibrium between hydrogen adsorption and desorption, a necessary condition for catalysis over HER.
Although these findings help our understanding of the interaction between transition metal carbides and hydrogen, future work is needed both experimentally, to prove those findings, and theoretically, to investigate the consequences of the observed trends, including the feasibility of CO 2 reduction over the surfaces of these carbides.

Conflicts of interest
There are no conflicts to declare.

Acknowledgements
FS thanks the Cardiff University School of Chemistry for a fully funded PhD scholarship. This work was funded as part of an