Surface oxidation for enhancing the hydrogen evolution reaction of metal nitrides: A theoretical study on Vanadium Nitride

The rational design of novel catalysts for energy applications has been an area of immense interest over the past decades. Recent studies have demonstrated that slight surface oxidation of nitrides can boost its catalytic activity, especially towards hydrogen adsorption. In this study, we have considered VN as a simple model nitride to simulate H adsorption on pure and oxygen polluted surfaces. First, we have evaluated the physical characteristics and stability of various facets. Diversity in the VN catalytic activity can originate from a variety of surface characteristics. Our calculations reveal the (112) and (200) surfaces as the most stable ones, while the (101), (111), and (113) surfaces are roughly in the same range of energy. The facet-dependent activity of VN toward H adsorption is then carefully discussed. A low level of oxygen contamination appears to robust the VN hydrogen adsorption and beyond that, a thin oxide surface layer can act as an activation layer, playing a positive role in improving the catalytic performance. We anticipate this picture to be an important input for designing enhanced nitride-based catalysts with controlled oxidation of surfaces.


Introduction
For energy-related challenges, many materials-oriented activities are in progress. For instance, there is an urgent need for developing highly active catalytic materials for processes involved in energy conversion and storage. There is an ever-increasing search for advanced rechargeable batteries with increased energy density, lower cost, long cycle life, and improved efficiency. Likewise, it is critical to develop H 2 adsorbents to meet the high volumetric/gravimetric density requirements for storing hydrogen while simultaneously maintaining thermodynamic-reversibility of the adsorption, with fast kinetics [1]. Water splitting and CO 2 reduction with low overpotentials and fast kinetics are significant areas of investigation as well [2]. Likewise, there are several efforts towards developing materials for supercapacitors to complement and improve other energy storage systems [3].
In energy-related materials research, transition metal nitrides have attracted considerable attention in recent years due to their potential as catalysts [4]- [6], hydrogen absorbents [7], and super capacitance electrodes [8]. Nitrides are very attractive catalysts for some reactions, including water splitting and NH 3 synthesis, because their catalytic properties which are comparable to those of Ptbased metals [9].
In recent years vanadium nitride (VN) has emerged as an interesting catalyst for many energyrelated systems. Vanadium nitride nanowires have been used to make gas diffusion cathodes of Li-CO 2 batteries [10]. VN shows promise as an electrocatalyst for both oxygen (OER) and hydrogen (HER) evolution reactions in water splitting, owing to its high availability, good electrical conductivity, and corrosion resistance [11]. Besides, there are several reports on hydrogen storage using VN-based nanocrystals. The reports indicate that the VN structure and its particle size are the most important factor for improving the gravimetric, thermodynamic, and kinetic properties associated with hydrogen storage [7]. Interestingly VN shows remarkable characteristics like high specific capacitance and a large potential energy window, along with benefits such as a good electric conductivity and a good thermal stability relevant for use as electrodes in capacitors [12,13]. However, the catalytic activity, stability, and adsorptive characteristics of VN are affected by the stoichiometry and crystal structure near the surface [12,14,15]. HER, OER, and gas diffusion on cathodes are structure-sensitive reactions [16]. The effect of particle size on the efficiency of hydrogen adsorbing mediums is inevitable [17,18]. Most redox reactions take place on or near the electrode surface of the supercapacitors [19]. And finally, in all these systems, the adsorption of hydrogen on the surface would be one of the main leading reactions. Therefore, a molecular-level understanding of the surface and its hydrogenation reaction stays a highly important task for the development of the aforementioned systems.
The existence of oxygen is unavoidable in nitride structures. According to some reports, oxygen incorporation with nitrides could be a deactivation factor, affecting their catalytic performance [20]. Interestingly, numerous reports have also shown the positive role that oxygen could play. A controlled oxide layer on top of the vanadium nitride surface leads to a superior capacity in electrochemical capacitor application of VN, without any loss in the overall electrical conductivity [19]. Furthermore, our previous work shows the excellent OER activity of the Nickel molybdenum nitride could be originated from a surface oxide-rich activation layer (SOAL) [21]. Recently, it is reported that VN capacitance could be even more increased by the presence of a single form of vanadium oxides, especially VO 2 , on the surface [15].
Experimental analysis of the surfaces is challenging and measuring the interaction strength between hydrogen and catalyst surfaces are difficult. Significant developments have been made in the use of density functional theory (DFT) for accurate simulation of the surface catalytic reactions [22]. DFT offers a means to investigate the surfaces, simulate reactions, and determine the magnitude of the interactions among the reaction intermediates and the catalyst surface. The other well-known descriptors of catalytic performance, also accessible by DFT methods, are the d-band center and the geometric structure of catalyst surfaces. The d-band model is based on the electronic structure of the material and is defined as the average energy of d states when projecting the molecular orbitals onto the metal atom's d-state at the surface [23,24]. An effective strategy for enhancing the catalytic activity of a surface toward HER is tailoring its d-band center, adjusting the free energy of H adsorption (ΔG H* ) close to the optimal value for the HER (i.e., zero free energy of adsorption) [25]. Manipulating the d-band center position has been successfully applied for improving OER and HER performance of metal nitrides previously [26]. Achieving further development in the catalytic performance of nitrides could be aided through computational insights.
To the best of our knowledge, there is no prior study on VN surface structures despite the interest in VN for this wide range of catalytic applications. Investigations on hydrogen adsorption of VN and the effects of oxidation on it are also very limited, especially from a theoretical perspective [15]. The current study aims to provide a detailed and comprehensive DFT investigation of VN surfaces, hydrogen adsorption on pure and oxidized surfaces, along with offering in-depth information about the catalytic activity of this nitride.
In this study, we first consider the stability of VN facets with both low and high Miller indices by calculating the surface energies and investigating the surface physical properties. Then the H adsorption is studied on clean facets and also the surfaces with different rates of oxygen contamination. Finally, a single layer of oxygen is assumed on top of the VN surfaces and various coverages of H will be discussed below.

Computational method
We have used density functional theory calculations in the generalized gradient approximation [27], as implemented in the Quantum Espresso package [28]. The projector augmented wave (PAW) [29] and optimized norm-conserving Vanderbilt (ONCVPSP) [30] pseudopotentials are utilized to describe the electron−ion interactions. The minimum energy cutoff is 45 Ry for the wavefunctions and 400 Ry for the charge density. A 6×6×1 Monkhorst-Pack k-point mesh was used for sampling the Brillouin zone.
VN in bulk structure has a halite, rock salt crystal structure in the cubic Fm-3m space group. The calculated lattice parameter is a= 4.28 Å for a fully relaxed structure. Each V atom has six equivalent bonds with the adjacent N atoms, where the bond length is 2.06 Å.
We have conducted DFT calculations for eleven cleavages on eight VN facets to evaluate the stability of the crystal facets. The surface calculations are accomplished by a slab model, while the primary unrelaxed surfaces are set up using the relaxed bulk structure. The basis of the bulk unit cell and the atomic positions are transformed such that the interested Miller index is oriented parallel to the (001) plane [31]. The slabs should be thick enough for convergence calculations, so we have considered the twelve-layer-thickness slabs. A vacuum of 15 Å along the c direction is sufficient for separating the periodic images. The lattice parameters of the slabs are fixed at the VN bulk value, while the atoms on the upper half of the slabs are allowed to relax until the residual force on each atom is smaller than 0.001 eV/Å. The lower half of the slabs are frozen to simulate the bulk structure. Since we have kept atoms fixed on one side, the slab has two inequivalent surfaces. Therefore, the unrelaxed surface energy, which is the energy of the surface before any surface optimization, is also considered for finding the final surface energy of the relaxed surface [32]. The surface energy (γ) is calculated by the following formula, Where the second term denotes the surface energy before any surface optimization. E b , E u , and E r resemble the total energy of the primitive cell, the unrelaxed slab, and the relaxed slab with one side fixed in the optimized bulk positions, respectively. The number of atoms in the slab in comparison with the ones in the bulk is shown by n. A is the surface area of one side of this slab.
For simulating the H adsorption of the surfaces, we have adopted a procedure similar to the HER simulations [33]. HER occurs in two steps, as represented in equation (2) where the (*) symbol represents a free adsorption site and H* denotes an adsorbed hydrogen atom on a facial adsorption site. Different H adsorption sites are considered where the initial distance of H to the surface is 3 Å. The Gibbs free energy change ΔG for the hydrogen adsorption is considered as an effective and powerful gauge for indicating the activity of HER catalyst [23], [33]. H* is the intermediate agent for both proposed mechanisms of HER in equation (3), i.e. the Volmer-Heyrovsky and the Volmer-Tafel mechanisms [34] and ΔG can be widely used for describing the HER electrocatalytic activity [24], [35].
The Gibbs free energy change is calculated as is the difference between the calculated reactants and product energies: ∆E ∆ = ( ( 2 ), (5) where and represent the surface total energies with and without the hydrogen adsorption, respectively. The number of adsorbed H atoms are shown by n, and E(H 2 ) denotes the total energy of a hydrogen molecule in the gas phase. denotes the zero-point ZPE energy, is the variation between the ZPE of the adsorbed hydrogen (H*), and H in the gas ∆E ZPE phase: T is the room temperature and is the entropy changes: ∆S where stands for the entropy of H 2 gas at normal conditions. According to the Nørskov et al.

S(H 2 )
calculations, the vibrational entropy in the adsorbed states is negligible, and is considered as T∆S 1 2 equal to -0.204 eV at room temperature [33]. If the ΔG H* has a negative value (equation TS (H 2 ), (4)), H can be easily adsorbed on the surface, while its desorption would be difficult. Therefore, the ideal HER electrocatalytic activity will belong to the material whose ΔG H* value is zero; while the energy range between -0.2 to 0.2 eV is normally considered as the criteria for showing active HER catalyst for a sample.
Oxygen doping is considered by replacing a nitrogen atom with oxygen for simulating a minor Odoping, and two out of the eight N atoms of the slab's uppermost level for simulating an oxynitride surface with 12.5% oxygen on top. Then, an oxide layer is put on the VN surface and all N atoms are replaced with O atoms. This results in a VO/VN configuration. The presence of impurities leads to a surface non-uniformity, which increases the diversity of possible hydrogen adsorption sites from a theoretical point of view. Lastly, the coverage of hydrogen atoms on the surface was raised from 1/8ML to 1/4ML and 3/8ML on pure and VO/VN surfaces.

Surface investigation
Different synthesis methods would result in different multi-faceted samples. Finding the reactive facets and explaining their role in promoting the catalytic activity of a material at a molecular level is highly valuable. This information not only provides a guideline for uncovering the catalytic mechanisms behind, but also can lead to the design of new catalysts.
The most basic step is evaluating surface stability. Eight low and high-index VN surfaces have been selected from both types of flat and step/edge facets. We have considered four previously synthesized VN surfaces, i.e. (200), (111), (220), and (311) [7,15], while (001), (112), (113), and (101) surfaces will be discussed as a predictive model. Fig. S1, shows the lateral view of all considered slabs, before and after the structural optimization, and Table S1 reports the surface energies. Lower surface energy indicates a more stable surface. The stability sequence is (200)> (112) > (101), followed by (111), (311), and (113) in close stability with each other. The (220) and (001) facets have the highest surface energies and accordingly are the most unstable considered surfaces. A brief discussion on the physical properties of the surfaces is given below: The VN (001) and (220) are flat facets with rather smooth structures. A 4-fold atomic arrangement exposes on top of the (001), similar to all fcc structures. Both types of V and N atoms are present on the surfaces, exposed to vacuum. After the surface optimization, there would be a slight decrease in some V-N bonding lengths, which would be compensated by increasing the length of adjacent bonds. In general, the changes due to the surface effects are negligible on both facets and as was mentioned above, these facets have the highest surface energy among all considered facets.
The (311) surface is observed on the VN sample with an oxide surface [15]. Due to the geometric similarity of (113) with (311) surfaces, the geometric structure and stability of (113) are also examined in this study. At the atomic scale, two different cleavages can be considered on top of these step facets: the slabs with a V-rich layer on the topsides (referred to as the (113) (Table S1) show Vterminated cleavages are more stable on both (113) and (311) with almost the same surface energies. To reduce the effect of surface dangling bonds, the topmost V atoms move toward the bulk level, ending to relatively smooth facets.
The VN (111) crystal plane is also experimentally observed during the HER [16,35], and OER [11,16] studies. Besides, VN (111) has been considered as a substrate nanosheet for improving the OER performance, providing very desirable results [11,36]. In all fcc crystals, the (111) layers are closepacked. On VN (111), consisting of two interpenetrating fcc structures, all surface atoms have a relatively high coordination number. In accordance with (113) Table S1, respectively. After the geometric relaxation, atoms in the direct vicinity of the vacuum move downward, toward the bulk layers. This leads to the shrinkage of V-N bond lengths of up to ~9%, to neutralize the unphysical forces of the surface. The stability of this surface is close to the ones of (113) and (311). However, there is a small energy difference between the (111)-I and (111)-II cleavages, so there is a possibility for the pure nitrogen surface to be the uppermost layer.
In the following, we will show that N sites are not catalytically active sites toward hydrogen, which can be a major drawback for the catalytic activity of (111) surface.
The VN (112) surface is a completely open surface. Different colors are used in Fig. S2(a) to display the different heights of atoms relative to the vacuum (Z-axis). The atoms on the second and thirds layers are also visible through the holes of the first layer, although these atoms would not be directly accessible for an external gas molecule. Both N and V are in the vicinity of the vacuum and the difference in height of atoms in the vicinity of the vacuum are up to 4Å before the surface relaxation. Significant changes occur on the surface due to the surface effects, so that the height and bonding lengths of the uppermost atoms are substantially decreased. This is the most energetically favorable surface, which will be discussed more in the next sections of the paper.
(200) is the second most stable surface according to our calculations ( Fig. S1 and Table S1), which is in agreement with the experimental reports, synthesizing the uniform VN (200) facets [15,36]. Surface effects cause the atoms of the first layer to move downward. This increases the dependency of the V atoms to the lattice, while giving more freedom to the topmost N atoms, in a way that the final relaxed surface is slightly rippled (Fig. S1(k)). The surface reactivity can be strongly dependent on the activity of the terminating layer.
In the following, we will confine our further studies to the VN (112) and (200) surfaces, which have the lowest surface energies. This selection is not only based on these energy considerations, but also on the significant geometric changes at the surfaces, the existence of both types of N and V atoms at the surfaces, and other reasons which will be discussed below.

Gibbs free energy and H adsorption profile on clean VN surfaces
A rational design of novel catalyst materials for energy applications is essential. However, this work is challenging because of the complexity of the mechanisms involved. A large number of experimental and theoretical guidelines have already been established to find optimal catalysts [24]. One of the most successful theories is known as the Sabatier principle [38], which suggests that the binding strength between a catalyst and adsorbates is the leading descriptor of the catalytic activity. In this framework; a moderate interaction leads to the best performance, while too weak or too strong interactions would limit the adsorption or desorption rate of the catalyst surface, respectively.
Both VN (112) and (200) surfaces offer a variety of possible sites for adsorbing gas molecules. These sites can be considered in two general categories: i) On-top sites, where H atoms are initially considered on top of the metal or nitrogen atoms. ii) The hollow-sites, where H is considered above the vacant space among four atoms. We are not aware of any specific report showing the catalytic activity of VN as a function of adsorption sites.
We have considered two dissimilar positions of nitrogen for H adsorption (with two and three dangling bonds), one position on top of V, and one bridge site on the VN (112) surface. Fig. S2 shows the schematic view of the surface, with the H* adsorption sites marked as site L1-4. The L2 and L3 sites, locally shown in Fig. S2(b), are the most active sites according to the Gibbs free energy results (Fig. S2(f)). Here, the H atoms are absorbed by the V atom with a moderate V-H bond length. On the other hand, the H atom on the N site (noted as L1) is strongly adsorbed and a short H-N bonding, with an approximate length of 1.02Å, is created on top of the surface (Fig.  S2(a)). As expected, the most negative free energy, and accordingly the lowest H evolution activity, is related to this site (Fig. S2(f)).
Similarly, on the (200) surface there are three main possibilities for H adsorption, noted as A1, A2, and A3, referring to H atoms on top of the V, N, and hollow sites, respectively (Fig. S3). The bridge site (A3) is the most stable adsorption site, which is followed closely by A1. Their related free energies are clearly lower than the corresponding sites on (112) surface. On the other hand, a short strong H-N bonding on A2 indicates that the site won't be active toward HER. Both A1 and A3 sites have reached a similar final structure after surface rearrangements ( Fig. S3 (a) and(c)). The distance of the hydrogen-bonded V atom is considerably increased with respect to the other V atoms and also the neighboring N atoms. The lengths of H-V bonds are around ~1.67Å on the most active sites of the clean (200) surface.
Therefore, the best adsorption sites of the clean surfaces are the topmost vanadium atoms. However, the V-H bonds are relatively strong and the corresponding free energies toward HER activity are quite negative (-0.52 eV for A3 and -0.55eV for L2 according to our results). In the following, we examine whether oxygen impurities can yield a positive effect on the catalytic performance of the surfaces.

Oxygen contamination of VN surface
The H adsorption of VN is studied while the surfaces contain different amounts of oxygen impurities. Four contamination ratios have been studied, where 12.5%, 25%, 50%, and 100% of the N atoms on the topmost layer have been replaced by oxygen. Fig.1-2, Fig. S4-S5, and Table 1 show the corresponding relaxed structures and H adsorption free energies.
We will first discuss the O-doped VN (112) Fig. S4(a)). Interestingly, all three latter situations (i.e., M5, M6, and M7) lead to the same final surface configuration (Fig. S4(c)). The adsorbed H atom forms two bonds with the lengths of 1.93-1.95 Å with V atom, while is placed 3.15 Å away from oxygen. As expected, the Gibbs free energies of M5, M6, and M7 sites are energetically very close ( Fig.  S4(d) and Table 1 Fig. 1(j) by B1-6. The local configurations around H adsorption sites are represented in Fig. 1(a)-(g). On four sites, i.e., B1, B3, B5, and B6, adsorbed hydrogen atoms are connected to vanadium. The V atom in the immediate neighborhood of O is the most active site due to the energy results reported in Table 1. and Fig.  1(k). In comparison with the clean surface, the Gibbs free energies are more favorable for HER. Moreover, on all the aforementioned sites, the V-H bond lengths are a bit longer than the clean surface. On the other hand, H is strongly absorbed by the oxygen atom on B2 site ( Fig. 1(b)), constructing an OH unit with a length of 0.978 Å. Later, we will discuss about this unit more. Both H-N bonds ( Fig. 1(d) and (g)) are also short and strong, therefore, oxygen impurity does not contribute to the better catalytic activity of these sites.
By increasing the surface contamination rate, we find that V atoms in the direct vicinity of O atoms still would be the most active sites of the surface. In Fig. S5, 25% of the N atoms on the outermost layer are replaced with oxygen. Not only the H adsorption energies of C1, C2, C4, and C5 are less than that of the clean (200) surface, but also the bonding length between V and H are significantly increased (Fig S6(b)-(h)). Interestingly, there is no OH unit formed on the C1 site. However, despite the longer H-V bond lengths on this site, its related free energy (~0.47 eV) is more negative in comparison with C2 and C4 sites, with ΔG H* around -0.42eV and 0.47eV, respectively. It is due to the existence of two H-V bonds at this site, which will increase the dependency of H on the surface and thereby reduce the catalytic activity.
The efficiency of a catalyst can be better understood in terms of its electrical properties. According to the density of state (DOS) plots for the pure VN surfaces (Fig. S7(a), Fig. S8(a)), the interaction between V and N is expectedly strong. The transferred charge between vanadium and nitrogen results in the DOS being similar to that of the Pt-group metals near the Fermi level. This similarity shows why VN has been considered as a potential candidate for catalytic activity toward evolution reactions. In doped structures, oxygen will influence on the vanadium and nitrogen electronic states ( Fig. S7(b)-(c), Fig. S8(b)-(d)). The main peak associated with oxygen is in the energy range between -9 and -4 eV, where the creation of V-O and also N-O bonds are clear. By increasing the rate of oxygen impurities, this peak shifts toward the Fermi energy. The other important O-p DOS contribution is around the Fermi level, with -4 to 5 eV energy. The V-d and O-p orbitals are degenerate in this area, resulting in both bonding and antibonding states. However, the O-p contribution is so less than that of V-d states, so the d-character remains near the Fermi level. By increasing the rate of surface oxidation, the total density of states is even slightly increased near the Fermi level ( Fig. 3(c), which promises a better catalytic activity.
Tailoring the d-band centers is another good way for qualifying nitrides as highly active catalysts for HER [26]. The d-band center is a physical parameter, strongly correlated to the adsorption energy of the adsorbate [33]. By increasing the oxygen impurity rate, the VN d-band center shifts toward the Fermi level according to Table S2. The higher value of the d-band energy represents stronger adsorption [25]. Thus, this shift implies a stronger H adsorption to the VN surface and hereby the lower value of the free energy (ΔG H* ).
In summary, the surface energy of VN (112) and (200) facets are in the same range, but both clean and O-doped VN (200) surfaces show a better catalytic activity toward HER according to our results. This result could be so helpful in practice, providing a direction for laboratory synthesis. Besides, oxygen contaminations of the facets cause the enhancement of the catalytic activity of the polluted surfaces. The active sites are V atoms in the direct neighbor of oxygen. This easy surface engineering would result in modulation of the d-band center and could be a good strategy to robust the hydrogen evolution. The aforementioned reasons persuade us to expand the oxygen coverage, placing one oxide layer on VN (200) surface, and investigate the catalytic activity of this framework.

VO/VN configuration and catalytic activity on different H coverages.
Various vanadium oxides have also been reported favorable for HER, but their low inherent conductivity and insufficient cycling performance are the main issues for the practical applications [39]- [41]. By considering VO/VN configurations, the good electrical conductivity of the VN in combination with the vast variety of active hydrogenation states could boost the nitride's catalytic activity [19]. Besides, we will benefit from a protective oxide layer for enhancing the stability of the material.
After replacing all uppermost N atoms with oxygen, we obtain a single vanadium oxide layer on the VN (200) surface. DFT structural relaxation directs to V-O bonds with average lengths of 2.08 Å on the VO/VN surface. The possible H* adsorption sites (D1 to D7) and its lateral view are shown in Fig. 2(h). The oxygen atoms are in three non-uniform surface positions, which are marked as D1, D3, and D5 in Fig. 2. The D1 site has the lowest position among than all other topmost atoms. It moves away from adjacent atoms after structural relaxations and the only remain concrete bond is with a V atom on the following layer. The free energy values are reported in Table 1 and a schematic comparison between the best adsorption sites is shown in Fig. 3(d). The results clearly demonstrate the reduction of energy amount after oxygen doping of the surfaces, especially when there is a thin layer of oxide on top of the surface. The energy values for VO/VN surface follow this order: D1< D2 ≈ D3 < D4 < D6 < D5 ≈ D7, where the best energy values are -0.32 eV for D1, -0.36 eV and -0.38 eV for D2 and D3, respectively. The value of ΔG H* for these sites shows a considerable reduction (in terms of energy amount) compared to the oxygen-free surfaces (with the value of -0.52eV for A3 and -0.55eV for L2 on (200) and (112) surfaces, respectively). On the contrary, the difference between the ΔG H* of D4, D5, D6, and D7 sites with that of the to the clean (200) surface is vital (Table 1), but the V-H bonds on these sites are stretched out up to 13% ( Fig.  2(b)-(d)), which shows the sites could be also more active in comparison with the clean surface.
More interestingly, the H atom is easily adsorbed by adjacent V atoms of D1 and D5 sites ( Fig.  2(a), (c)), while the H atom on top of the D3 site is absorbed by oxygen and forms an OH unit. Oxygen in the D3 has longer V-O bonds and less dependency on the surface, which could easily bound to H and freely desorbed from the surface. The Metal-OH bond strength is usually a theoretical descriptor for OER, meaning that a balanced bond strength between the metal atom and OH unit leads to optimal activity [41,42]. Given that the OH unit has left the surface freely ( Fig.  2(g)) and there is no strong bonding with the V atom, it does not seem like an active OER site.
Although we do not calculate the binding energies of OH* on the surface and do not intend to simulate OER, in the following paragraph we have demonstrated that higher H coverage will result in the creation of an OH unit with the appropriate dependency to the surface on site D5, which may show good OER activity of VO/VN surface.
Finally, we have increased the coverage of hydrogen from 1/8ML to 1/4ML and 3/8ML on both VN (200) and VO/VN surfaces. We have considered different adsorbing sites of the second hydrogen on the clean surface as shown in Fig. S9(a). We have supposed that the first H atom is adsorbed by one of the top-most V atoms due to our previous discussions. The second H will be adsorbed by another V atom ( Fig. S9(b), (d)), or by an N atom (Fig. S9(c), (e)), or forms a hydrogen molecule and leaves the surface (Fig. S9(f)). As this process continues and the surface hydrogen coverage increases, the V-H bonds are easily broken and the hydrogen molecules freely desorb from these sites ( Fig. S9(g), (h)). However, the N-H bonds are short and strong, and the bond breakage occurs only in high H coverage. Therefore, it is unlikely to have an H occupied surface due to the high activity of the V sites, although the strong N-H bonding will decrease the surface HER activity.
On VO/VN surface, the second H is also adsorbed freely by another V atom, and two equal V-H bonds are formed with lengths of 1.677 Å. After the structural relaxation, a small hole is created on the surface, which is clear by a non-bonding area in ELF plots, displayed in dark blue in Fig.  4(c), (d). Considering the small size of the H atom, the presence of this hole can create an easy path for hydrogen to penetrate to the underlying layers and enhance the HER activity.
By adding the third H near the adsorbed hydrogen atoms, as expected, the hydrogen molecule detached from the surface (Fig. S10). It indicates that the sites are fully active for HER. However, in another interesting case, when we have placed the third hydrogen near O atom on the D5 site, an OH unit is formed on the surface. The unit is connected to three topmost V atoms. The V-O bonding lengths are 2.0-2.2Å, and the bonds are neither too fragile nor too strong following the ELF results ( Fig. 4(d)

Conclusion
Nitrides have been investigated as potential replacements for noble metal electrocatalysts in a wide range of applications. One of the most effective factors on their catalytic activity is the surface capability for binding or reacting with adsorbates.  Fig. 1), and double O-doped (C sites are shown in Fig. S5) VN (200). D sites are used to represent the sites on VO/VN (200) surface (Fig. 3), where L and M sites are on VN (112) surfaces shown in Fig. S2 and S4 Table 1. An asterisk is placed on some sites to indicate that these atoms are linked to hydrogen.