Is the doped MoS₂ basal plane an efficient hydrogen evolution catalyst? Calculations of voltage-dependent activation energy

Abstract

Transition metal dichalcogenides are cheap and earth-abundant candidates for the replacement of precious metals as catalyst materials. Experimental measurements of the hydrogen evolution reaction (HER), for example, have demonstrated significant electrocatalytic activity of MoS₂ but there is large variation depending on the preparation method. In order to gain information about the mechanism and active sites for the HER, we have carried out calculations of the reaction and activation energy for HER at the transition metal doped basal plane of MoS₂ under electrochemical conditions, i.e. including applied electrode potential and solvent effects. The calculations are based on identifying the relevant saddle points on the energy surface obtained from density functional theory within the generalized gradient approximation, and the information on energetics is used to construct voltage-dependent volcano plots. Doping with 3d-metal atoms as well as Pt is found to enhance hydrogen adsorption onto the basal plane by introducing electronic states within the band gap, and in some cases (Co, Ni, Cu, Pt) significant local symmetry breaking. The Volmer–Heyrovsky mechanism is found to be most likely and the associated energetics show considerable dopant and voltage-dependence. While the binding free energy of hydrogen can be tuned to be seemingly favorable for HER, the calculated activation energy turns out to be significant, at least 0.7 eV at a voltage of −0.5 V vs. SHE, indicating low catalytic activity of the doped basal plane. This suggests that other sites are responsible for the experimental activity, possibly edges or basal plane defects.

---

1 Introduction

Hydrogen is considered as one of the most promising means of future storage of renewable energy from intermittent sources.¹ Gaseous hydrogen can be produced through electrolysis of water in a potentially cheap and sustainable way of converting electrical energy from renewable sources into chemical energy. The efficiency and cost of this process largely depends on the catalyst material. Currently, the process relies heavily on the high activity of platinum-group metals (PGMs), and particularly platinum itself.² The limited availability of these precious metals in Earth's crust and socioeconomic issues in the mining countries pose problems for long-term sustainability and calls for the development of new catalyst materials as well as better understanding of the fundamental mechanisms of the hydrogen evolution reaction (HER) and electrocatalysis in general. Among the emerging candidates in the search of sustainable replacements, various metal alloys and transition metal compounds have been found,^3–6 with recent emphasis on engineering of two-dimensional materials.^7,8 Transition metal dichalcogenides have been shown to exhibit promising properties for HER, and MoS₂ has, in particular, been widely investigated both experimentally and theoretically.^9–13

An important step in characterizing and designing new electrocatalysts is to identify the atomic sites that exhibit high activity. While the edge sites of pristine 2H-MoS₂ have been shown to be catalytically active, the basal plane is inert in its pure form.^9,14,15 MoS₂ has been synthesized in a wide range of morphologies with a focus on increasing the abundance of edge sites.^12,16–18 Further, the basal plane can be activated by introducing defects such as S-vacancies,¹⁹ phase boundaries²⁰ and metal²¹ or non-metal²² impurities. As such, the full picture of MoS₂ activity is quite complex and a basic understanding of the features in atomic and electronic structure that enhance HER is needed to optimize performance.

Experiments have indicated that transition metal doping enhances the overall activity,^23–26 but detrimental effects and conflicting results have also been reported,^27,28 illustrating that the result strongly depends on the system specifics, notably morphology, the nature and level of doping, and the experimental techniques applied. This indicates that the manifested activity relies on the interplay between several factors and that theoretical studies can therefore be helpful for identifying the contributing ones.

In the present study, we investigate the effect of transition metal doping on the electrocatalytic activity of the MoS₂ basal plane. Theoretical studies have shown that hydrogen adsorption onto the basal plane is enhanced by transition metal doping,²⁹ and in the following we will assess whether this corresponds to higher activity (reduced reaction barriers) under electrochemical conditions for the whole sequence of 3d-metals as well as platinum. First we investigate the doped material itself and hydrogen adsorption in the gas phase. Then we move on to model the electrochemical reactions involved in hydrogen evolution. These are more challenging than calculations of gas/surface reactions because: (i) the reaction occurs in the presence of an electrolyte and (ii) the reaction occurs at a fixed electrode potential. These challenges must be overcome to provide a realistic model for electrochemical reactions, as further discussed in the following section.

2 Methods

2.1 Calculations of activation energy

All results were obtained from spin-polarized density–functional theory (DFT) calculations within the generalized gradient approximation (GGA). The revised Perdew–Burke–Ernzerhof (rPBE) exchange–correlation functional by Hammer et al.³⁰ was used because of its improved results for adsorption energy. In addition, van der Waals interactions were accounted for by the zero-damping D3 parameters³¹ as this provides an improved description of the interlayer distance of multilayer MoS₂. A cutoff of 400 eV was used for the plane wave kinetic energy in the representation of the valence electrons, and the projector augmented wave (PAW)³² approach was used to represent the effect of inner electrons. For all transition metals, the outermost s- and d-electrons were treated as valence electrons. For O and S, the 2s²2p⁴ and 3s²3p⁴ electrons were treated as valence, respectively. Test calculations including also 3s- and 3p-electrons for the early 3d-metals did not indicate any discrepancy in adsorption energies nor in the local electronic structure at the adsorption sites. As GGA tends to excessively delocalize the wave function due to self-interaction, the description of d-states was compared to that of the Hubbard U approach³³ and a hybrid functional (PBE0^34,35), see Fig. S4 and S5 (ESI†). The Vienna ab initio simulation package (VASP) was used in the DFT simulations.³⁶

The geometry optimization of bulk MoS₂ was performed using the primitive unit cell of 2H-MoS₂, with the Brillouin zone sampled by a Monkhorst–Pack (MP) grid of dimension 9 × 9 × 5. For calculations on doped mono- and bilayers we consider 5 × 5 × 1 supercells (75 atoms), with 3 × 3 × 1 MP grids. For accurate density of states (DOS) calculations, 11 × 11 × 1 MP grids were used. A vacuum layer of ca. 14 Å was introduced to decouple the periodic images of the slab.

The activation energy for the various elementary steps in the electrochemical reaction was calculated by first finding an approximate minimum energy path for the transition using the climbing-image nudged elastic band (CI-NEB)³⁷ method, followed by calculations with a tighter convergence as the saddle point on the energy surface corresponding to a given applied voltage is found using the minimum mode following (MMF)³⁸ method. The tolerance for force convergence in saddle point searches was set at 0.05 eV Å⁻¹ while the tolerance in minimization calculations was 0.02 eV Å⁻¹.

2.2 Reaction mechanism

In acidic conditions, the hydrogen evolution reaction involves the adsorption of H⁺ from solution onto the catalyst surface and subsequent desorption of gaseous H₂. This process can be described in terms of three steps:

H⁺ + e⁻ → H* (1)

H⁺ + H* + e⁻ → H₂ (2)

2H* → H₂ (3)

where H* indicates hydrogen bound to a surface site. Adsorption occurs through the Volmer (1) mechanism, and desorption through either the Heyrovsky (2) or Tafel (3) mechanisms. Whether the evolution proceeds through the Tafel or Heyrovsky mechanism (or a combination of both) depends on the kinetics of these in the given system.

The most widely used descriptor for the HER efficiency of a material is the free energy of hydrogen adsorption ΔG_H on its surface. In accordance with the Sabatier principle, a value of ΔG_H ≈ 0 has been shown to correlate with high exchange currents.^39,40 In the gas phase, we define the n-th hydrogen adsorption energy as

Furthermore, the Gibbs free energy is then given in terms of this energy as

ΔG_H = ΔE_H + ΔE_ZPE − TΔS_H. (5)

We can approximate the entropic term as where we neglect the configurational and vibrational entropy of the adsorbed state. At 298 K, the entropy contribution is about 20 meV. The zero-point energy E_ZPE is calculated individually for each system, and compared to the reference value of H₂ (vib. frequency of ≈4400 cm⁻¹), so that . We assume that the zero-point energy does not vary considerably for different hydrogen coverages.

Under electrochemical conditions, the chemical potential of a H⁺-e⁻ pair in solution is given in terms of the electrode potential U (V vs. SHE) and pH⁴¹ as

where k_B is the Boltzmann constant, T is the temperature and a_H+ is the activity of protons which is related to pH as pH = −log a_H⁺. Since we consider acidic solutions (pH → 0 corresponding to 1 M [H⁺]), the final term is negligible, and the chemical potential is essentially linearly modulated by the electrode potential.

2.3 Solvent model

One of the challenges of modelling electrochemical reactions is to ensure an accurate description of the solvent. Explicit description of the solvent is computationally intensive and the solvent atoms introduce a large number of degrees of freedom which complicates the process of finding energy minima and saddle-point structures. In this work, we employ an implicit solvent model through the implementation in VASPsol,^42,43 where the solvent (here water) is treated as a polarizable continuum. In this theory, the dielectric permittivity is spatially modulated between the extreme values of 1 (in vacuum) and ε_b (in bulk solvent) by a shape function ζ(r) as ε_r(r) = 1 + (ε_b − 1)ζ(r). The modulation defines regions of solute and solvent, depending on the local electron density n(r) in terms of the complementary error function. The free energy is minimized by equating the variation with respect to both the electron density and the electrostatic potential to zero. The former leads to additional terms in the local potential of the Kohn–Sham equations, as described in detail in ref. 43. From the latter we obtain the generalized Poisson–Boltzmann equation, which describes the distribution of counter-ions in solution.

The implicit solvent model allows us to largely omit explicit H₂O in the calculations, but to model a realistic stability and reactivity of the solvated H⁺, some water molecules are still needed. An often used model is the hexagonal ice bilayer. In most systems, however, the periodicity of the system of interest is not commensurate with the typical hexagonal water structure. We therefore employ a cluster model based on the Eigen cation (H₉O₄⁺), where a proton is effectively shared between four water molecules. For such a water model it follows that, since the ion is not externally restricted by hydrogen bonds, we expect to obtain a lower estimate of the activation energy compared to constrained structures such as hexagonal ice.

Test calculations were performed to assess the size effect of the water cluster. Calculations with Hydronium (H₃O⁺) and Zundel (H₅O₂⁺) ions showed respective discrepancies of at most ∼0.24 eV and ∼0.04 eV in the (Volmer/Heyrovsky) activation and reaction energies, with respect to the energies obtained with the Eigen cation (see Fig. S1, ESI†). This indicates reasonable convergence with increasing cluster size. We point out that calculations using only a single hydronium ion in the absence of an implicit solvent yield very different results for both the Volmer and Heyrovsky case. The reaction energies were decreased by roughly 1.1 eV and 2.7 eV, respectively, meaning that a single hydronium ion is far too unstable in the gas phase to provide reliable results. For the Eigen cation, the effect is smaller but still significant; the reaction energies are reduced by roughly 0.3 eV and 0.9 eV, respectively, in the absence of implicit solvent. Thus the implicit solvent is essential for the description of protonated water in combination with a cluster model.

2.4 Applied voltage

A finite simulation cell poses another challenge in determining the energetics of electrochemical reactions, as an electron transfer will lead to a substantial change in the electrostatic surface potential, and therefore to a capacitive contribution to the reaction energy. This contribution is inversely proportional to the lateral cell dimension, and the problem has previously been solved by extrapolating the energies to the limit of infinite lateral cell size.⁴⁴ This approach is quite expensive, and when dealing with constant dopant concentrations the cells rapidly become too large. Instead, we use a grand canonical approach (Fig. 1) where the electron number is allowed to vary to adjust the voltage.⁴⁵

Fig. 1 Example illustrating the methodology for converging reaction energy to a certain electrode potential, using parabolic fits. (a) The neutral (PZC) energy is compared to that at fixed electrode potential U = U^IS_PZC. The neutral calculations do not correspond to a specific potential, as it changes over the reaction with a magnitude depending on the size of the simulated system. The reaction is seen to become increasingly exothermic at lower potentials, with corresponding lowering of the required activation energy. For large negative potential, the transition state (TS) converges towards the initial state (IS). (b) Comparison of the obtained fits with the data points that are within a certain voltage threshold. The fit accurately reproduces the reaction and activation energy at a given potential. Data from the Heyrovsky reaction on Pt-doped MoS₂ with an initial hydrogen coverage of θ = 2/3.

Defining the electronic potential as μ = ε_F + eϕ_∞, where e is the electron charge, ε_F is the Fermi energy and ϕ_∞ is the electrostatic potential in the bulk solvent, the potential referred to that of the standard hydrogen electrode (SHE) is U = μ/e − ϕ_SHE. For the SHE we use the value ϕ_SHE = 4.43 V.⁴⁶ Further, we introduce the grand-canonical electronic energy Ω through the Legendre-transformation of the free energy f as Ω_U = f(n_e) + δn_eeU, where δn_e is the number of excess electrons in the cell.

Over a range in applied voltage that is not too large, the number of electrons is varied in increments starting from the neutral value, and the calculated energy of the system is then found to vary in a parabolic way with an extremum corresponding to the potential of zero charge (PZC). The curvature of the parabola represents the negative interfacial capacitance.⁴⁷ As illustrated in Fig. 1a, reaction and activation energies corresponding to a certain potential can then be obtained from the difference between fitted parabolas. The results shown in Fig. 1b coincide well with the alternative (and more computationally expensive, unless only a specific potential value is of interest) approach of converging each reaction geometry to the same potential, though extrapolating the fits outside the potential range of the available data cannot be considered reliable. Most systems display excellent parabolic behavior, but in some cases changing the electron number leads to changes in the geometry or crossing of small electronic gaps. This causes an abrupt change in the Fermi energy, and correspondingly in the potential.

3 Results

3.1 Doped MoS₂

The bulk 2H-MoS₂ phase consists of alternating layers, bound together by van der Waals forces, as shown in Fig. 2. By minimizing the energy with respect to the unit cell volume, we obtain an optimized structure with in-plane lattice parameter a = 3.18 Å, interlayer separation c = 6.20 Å, and layer height t = 3.13 Å. Both the hexagonal lattice parameters are within 1% of the experimental values (a = 3.15 Å and c = 6.15 Å, respectively⁴⁸). Our calculations show that pristine 2H-MoS₂ has band gaps of ∼0.9 eV (indirect) and ∼1.6 eV (direct) in multi- and mono-layer cases, respectively. This qualitatively corresponds to the experimental band structure, though the values are somewhat smaller due to the systematic underestimation of band gaps in the GGA approach. For comparison, the experimentally measured band gaps of bulk and monolayer 2H-MoS₂ are roughly 1.3 eV and 1.9 eV.^49,50

Fig. 2 Top: Structure of pristine 2H-MoS₂ with optimized lattice parameters. Bottom: Local geometries of MoS₂ doped with 3d-transition metals. Most dopants retain the six-fold symmetry, but notably Co, Ni, and Cu break this symmetry.

3d-Transition metal atoms are introduced as Mo-substitutional dopants to a 5 × 5 cell of the pristine monolayer, yielding a doping concentration of 4%. For Co, Cu and Ni, the local symmetry of the pristine geometry is broken, and the dopant atom binds to five surrounding S-atoms. This is due to excess d-electrons occupying antibonding orbitals.²⁹ For Sc, Ti, V, Cr, Mn, Fe and Zn, the symmetry is preserved (see Fig. 2). In the asymmetric case a sulfur atom has a broken or stretched bond, making it more exposed for adsorption (activated) with respect to the default coordination. We note that for Co, Ni, Cu, Pt, the dopant atom is also shifted down by 0.26–0.27 Å in the direction perpendicular to the plane (not visible in Fig. 2). For Zn there is a smaller shift of 0.15 Å. The remaining dopants are within 0.01 Å of the reference height.

Electronic states within the band gap are introduced as a result of the doping. Fig. 3 shows the density of states (DOS) of the various systems near the band gap. Apart from Cr, states are introduced in such a way that the effective band gap is significantly reduced. In the later 3d-transition metals, it appears that higher energy (occupied) S-p states are introduced, consistent with activation. The induced states are mainly localized at the sulfur atoms neighboring the dopant atom, and in the asymmetric cases the dislocated atom contributes the most. DOS and its projections are detailed further for surrounding S in Fig. S2 (ESI†). Since both the structural and electronic effects of the dopant are mainly local, larger distortions or reconstructions of the lattice are not expected.

Fig. 3 Projected density of states showing gap states emerging in 3d-metal doped MoS₂ monolayers. Gray lines indicate the Fermi level. Pristine MoS₂ in the bottom graph for reference.

3.2 Hydrogen adsorption

As mentioned above, the hydrogen adsorption free energy is a useful initial descriptor of the hydrogen evolution reaction, and therefore a natural starting point for investigation. We consider first the gas-phase situation, before moving on to solvated systems. Since hydrogen adsorption on the pristine basal plane is highly unfavorable, and due to the localization of the introduced S-p states, we expect a low surface coverage where only the dopant-induced favorable sites near the dopant atom (and neighboring sulfur atoms) are occupied. Fig. 4 shows the adsorption free energy for certain stable sites on Cu-doped MoS₂. Only the sites in immediate vicinity of the Cu-atom (A and B) are favorable, but since they can not be occupied simultaneously (see configuration D in Fig. S6, ESI†), we consider only the A-site for multiple H-adsorption in the following. This trend is similar in the remaining 3d metals, and the neighboring sites are lowest in energy also in the cases where adsorption is unfavorable. For details we refer to Fig. S3 (ESI†).

Fig. 4 Adsorption free energy of a single H atom coming from the gas phase and binding at various sites on Cu-doped MoS₂. The neighboring sulfur atoms are most notably activated. All sites display more favorable adsorption than the pristine basal plane.

The free energy is obtained from the electronic energy as explained in the previous section. The ZPE correction is calculated for a single H-adsorption for each system. Calculated E_ZPE for the dopant systems span a range of only 0.015 eV, and it can be approximated as a constant correction. In a similar manner, we assume that the vibrational frequencies do not change considerably as the coverage increases.

The free energy of adsorption for one to three H atoms in the simulation cell is given in Fig. 5a, where a coverage of θ = 1 refers to saturation of the three available A-sites. Consistent with the previous arguments of symmetry and electronic DOS, we observe that only the Co, Ni, Cu and Zn systems show energetically favorable adsorption of a single hydrogen atom, and higher coverages are even less favorable. Adsorption configurations that display moderate change in free energy are of interest for hydrogen evolution. This means that for systems Sc through Fe, only single H adsorption is relevant (note that Cr and Mn have high adsorption energies in any case). For systems Co through Zn, we consider an initial coverage of one or two H per dopant atom. For reference, the first adsorption energy was also calculated for the two-layer MoS₂ case of systems Fe through Zn, with only the top layer doped. Compared to the monolayer case, these values differ by less than 0.02 eV, suggesting that the hydrogen binding onto these dopant-activated sites is not sensitive to the presence of underlying layers. Monolayer calculations are thus appropriate for representing the general slab systems.

Fig. 5 (a) Incremental free energy of adsorption of multiple H-atoms for the doped and pristine MoS₂. (b) The promotion energy E_pdvs. free energy of adsorption for coverages up to the first unfavorable adsorption for each dopant.

A chemical picture of the adsorption mechanism on MoS₂ is that excess charge on the sulfur atom upon adsorption is partially distributed over the surrounding metal atoms.⁵¹ In this regard, the promotion energy is a useful descriptor for the adsorption energy: E_pd = E_d − E_p, where E_p is the p-orbital center of the S-atom in question, and E_d is the center of the unoccupied metal d-orbitals integrated from the Fermi level up to the point where one electron is added, i.e. the effective LUMO of the local system. Fig. 5b shows the correlation between the adsorption energy and E_pd. We note that the local geometry changes upon hydrogen adsorption, which leads to additional energy contributions and deviations from this descriptor. This is especially visible for the Fe-system, where the local six-fold symmetry is broken upon hydrogen adsorption. The other systems maintain their symmetry, but the distance between the dopant atom and the active sulfur changes. The promotion energy captures two important conditions that are necessary for favorable adsorption energy, namely that the sulfur p-orbital states must be high in energy (activated), and that the surrounding d-orbitals must have unoccupied states not far above the Fermi level. It follows that large band gaps are detrimental towards favorable adsorption in these systems.

3.3 Neutral cell hydrogen evolution

Next, the presence of the solvent is taken into account and the various steps of the hydrogen evolution are calculated where the systems are kept under neutral conditions, that is at the potential of zero charge (PZC). Afterwards the grand-canonical approach will be used to keep the electrode potential fixed during the reaction. We consider first the Volmer step, where a proton from the water cluster adsorbs on the MoS₂ surface. It follows a simple reaction path, depicted in Fig. 6 for the Ni-doped system at initial coverage θ = 1/3. As also demonstrated in Fig. 6, the relation between activation and reaction energies for the Volmer reaction on the doped systems follows the Brønsted–Evans–Polanyi principle, where the activation energy E is linearly proportional to the reaction energy. The Volmer step is highly unfavorable for all the early 3d-metals, even though Sc, Ti, and V have much lower adsorption energies in the gas phase. For Sc, Ti, V, Cr and Mo, there appears to be either no saddle point (transition state, TS) between the initial (IS) and final state (FS), or the TS is very close to the FS, as obtained with both MMF and CI-NEB searches. Thus, the activation energy tends towards the reaction energy as the reaction becomes more endothermic, and the Volmer reaction is effectively only uphill in energy in these systems. Importantly, the adsorption state is not kept (meta)stable by a reverse barrier and its lifetime is insignificant.

Fig. 6 Top: Reaction mechanism for the Volmer reaction step on Ni-doped MoS₂ at an intial coverage of θ = 1/3. Bottom: Calculated activation energy vs. reaction energy for the various dopants and coverages studied here at the potential of zero charge, and comparison with the Brønsted–Evans–Polanyi relation (dashed line).

The large difference between the gas phase adsorption energy and the Volmer reaction energy in solution obtained for some of the systems (in particular Sc, Ti, V) shows that the solvent plays an important role in the reaction energetics. That is, the water cluster interacts differently with the surface before and after its proton has been transferred. For example, for systems with favorable hydrogen adsorption, the cluster is weakly bonded to the adsorbed hydrogen through a hydrogen bond. Such a configuration is not stable for the systems with unfavorable hydrogen adsorption, as it would lead to the hydrogen atom re-entering the solution. The cluster is thus further from the surface in a metastable final state, not supporting the attractive interaction. The correspondence with gas phase adsorption becomes clear again if one considers the next step of the reaction. After the Volmer step, the water cluster is again supplied with a proton from the bulk acidic solvent. The difference in energy between this state and the initial Volmer state (correcting for the additional H atom), is what corresponds to the adsorption energy. However, what influences the Volmer barrier is the initial Volmer reaction energy, as is clear from Fig. 6 and the linear relation. This illustrates one aspect of the significance of the choice of solvent description. We note that there is a small energy barrier of roughly 0.15 eV associated with the proton transfer step from the bulk solution. This is an upper bound in the sense that proton diffusion towards a negative surface will be associated with a favorable free energy change given by the potential difference vs. the bulk solvent, thus also lowering the associated barrier.

In the Heyrovsky step, a proton in the water cluster binds to an adsorbed hydrogen atom, and an electron is transferred from the electrode to form an H₂ molecule that is released from the surface. The reaction path is shown in Fig. 7. As in the Volmer reaction, the activation energy follows a linear trend with respect to the reaction energy. In this case, the four last 3d-metals, as well as Sc and Ti, are seen to have highly unfavorable Heyrovsky reaction energies at low coverage. For the late 3d-metals this is expected from the Volmer energies, but Sc and Ti appear to have highly unfavorable energies for both reactions. This is related to the previous discussion on the distinction between the initial Volmer reaction and the following proton resupplying step. Once resupplied, Sc and Ti return to lower energies, while V does not. This is consistent with the obtained Heyrovsky energetics.

Fig. 7 Top: Reaction mechanism for the Heyrovsky reaction step on Ni-doped MoS₂ at an intial coverage of θ = 2/3. Bottom: Calculated activation energy vs. reaction energy for the various dopants and coverages studied here at the potential of zero charge, and comparison with the Brønsted–Evans–Polanyi relation (dashed line).

The Tafel step involves the desorption of two adsorbed H-atoms on the surface. The kinetics of this step therefore depend heavily on the geometries and relative energies of the adsorption sites, as well as the surface H-coverage. This step depends less on the presence of the solvent than the Volmer and Heyrovsky steps. Adsorption of several H atoms will naturally involve occupying the less favorable sites, as well as introducing H–H interactions. From the single-H adsorption energies, we expect that subsequent adsorption will still be confined around the dopant atom, although energetically less favorable due to H–H repulsion. As shown in Fig. S6 (ESI†), it is evident that the dopant-induced sites are still preferred despite repulsive H–H interaction. Referring to the labeling of the sites in Fig. 4, the optimal reaction path proceeds from a neighboring A–A configuration through the A–B and B–B intermediates. For details on the reaction path and resulting activation energies and scaling relation, we refer to Fig. S7 and S8 (ESI†).

Comparison of the activation energies of the Tafel step and the Heyrovsky step at the same hydrogen coverage (θ = 2/3) shows that HER is more likely to occur by a Volmer–Heyrovsky mechanism rather than a Volmer–Tafel mechanism, for all dopants except Zn. In these systems the barrier of the Tafel mechanism is significantly larger than for the Heyrovsky mechanism at the same initial coverage, but for Zn, the Heyrovsky and Tafel barriers are comparable, at 1.70 and 1.67 eV, respectively.

3.4 Constant electrode potential

So far, neutral systems at the PZC have been considered. We investigate next the energetics when the reaction occurs at a certain constant potential by varying the number of electrons as explained earlier. Fixing the electrode potential will lead to a certain correction to the PZC energies as demonstrated in Fig. 1a. The preferred reaction mechanism is not expected to change for most systems, but in general it will be a function of the applied potential. The Volmer–Heyrovsky process is expected to become more favorable the more negative the applied potential is, while the Tafel step is only weakly affected. Calculations for Fe through Zn confirm that the Tafel barrier potential dependence in general is weaker than that of the electron-transfer processes. For Zn_θ=2/3, the Tafel energy barrier is lower than the Heyrovsky barrier above −0.8 V, but such a high equilibrium coverage, θ = 2/3, is only reached at voltage below −0.5 V, so the Volmer–Tafel path will thus be (slightly) preferred in a range within this voltage window, but otherwise the Volmer–Heyrovsky path is preferred. For all other systems, the Volmer–Heyrovsky path remains more favorable at all potentials. With this in mind, we consider for simplicity the Volmer–Heyrovsky path of all systems in the following comparison.

Around 0 V vs. SHE, the adsorption energies are quite similar to those in the gas phase, and the configurations we noted earlier are still of interest. Therefore, we consider the initial hydrogen coverages of 0 (Sc, Ti, V, Cr, Mn, Fe, Mo), 1 (Co, Ni, Zn, Pt), and 2 (Cu) per dopant atom in equilibrium. For reference, the initial zero-coverage is also considered for all dopants.

Fig. 8a shows the potential dependence of the Volmer and Heyrovsky barriers for each investigated system in the range between 0 V and −1 V, evaluated at the equilibrium coverage at 0 V. Both barriers are lowered by the negative applied potential. Most notably, we observe that the Co- and Ni-doped systems with an initial hydrogen atom display significantly smaller Heyrovsky barriers than the rest of the systems with small Volmer barriers, especially around U = −0.5 V. Since balanced moderate barriers in general will lead to faster reaction kinetics than one small and one large barrier, these systems seem to have the most active basal planes for HER. Note again that to a first approximation, the adsorption energy is modified by eU (see eqn (6)) such that higher coverages will be favorable at large negative applied potentials. For U = 0 V and U = −0.5 V, endothermic adsorptions are included for all systems as seen in Fig. 8b, ensuring that the relevant coverages are considered. In the case of U = −1 V, the large potential could lead to a further increase of the H-coverage, which here would lead to qualitative increase (decrease) in the calculated Volmer (Heyrovsky) barrier. The pristine MoS₂ shows a particularly strong dependence on the applied voltage and surprisingly exhibits low barriers at large negative potential (U = −1 V).

Fig. 8 (a) Activation barriers for the Volmer and Heyrovsky steps at U = 0 V, U = −0.5 V and U = −1.0 V for the various dopants. The initial H-adatom coverage is taken to be the equilibrium coverage at U = 0 V. Dashed bars indicate the values for the case of θ = 0 for systems Co, Ni, Cu, Zn and Pt. Since the first adsorption is favorable in those systems, it is seen to be associated with a small Volmer barrier, and a correspondingly large Heyrovsky barrier. (b) Relation between the activity obtained from a simple kinetic model and ΔΩ_H for each studied dopant and coverage at U = 0 V and −0.5 V. Dashed lines show fits of the calculated points, while solid lines show the theoretical activity calculated from the relation between the Heyrovsky barrier Ω_h and hydrogen adsorption energy ΔΩ_H. The slope of this relation determines the resulting volcano shape.

A detailed evaluation of reaction kinetics would require calculations of several other parameters, such as pre-exponential factors. However, using the fact that all the systems studied are comparable, we can estimate the relative kinetics with a simple kinetic model. We characterize the turnover frequency (TOF) f by where is the probability of being in the state with i adsorbed hydrogen and Ω_h is the Heyrovsky barrier. is the partition function. At temperature T = 298 K, the resulting TOF assumes the characteristic volcanic shape with respect to ΔΩ_H, as seen in Fig. 8b. Points are calculated directly from the energetics for each system, while the solid lines represent the theoretical activity obtained by using the scaling relation between ΔΩ_H and Ω_h at the corresponding potential. This correlation is shown in Fig. S9 (ESI†), where the same outliers as in Fig. 8b are visible, notably Sc, Ti, Cr, Ni at U = 0 V and Sc, Ti, Ni at U = −0.5 V. These outliers lead to the different shapes between the two volcanoes. Peaks of the dashed volcanoes are shifted towards positive ΔΩ_H, indicating that in these systems the optimal condition is a slightly unfavorable adsorption rather than the perfectly neutral one, due to the competition between the Volmer and Heyrovsky barriers not being balanced at ΔΩ_H = 0 eV.

Further, we note that the transition state geometries change slightly as a function of the potential. For the Heyrovsky reaction, the distances r_H–H (between the two reacting hydrogen atoms), r_O–H (between oxygen atom and proton), and r_S–H (between sulfur atom and adsorbed hydrogen atom) largely define the reaction coordinate. As the potential is lowered, the magnitude of these values at the transition state tend to increase, decrease and decrease, respectively. In Fig. 9a, this is illustrated for all systems considered. Regardless of the dopant, a low Heyrovsky barrier is associated with a transition state which is geometrically more similar to the initial state than the final state, in that r_O–H is close to the H₃O⁺ bond length and r_H–H is far from the H₂ bond length. Thus, knowledge of the arrangement of the water molecules at the transition state is largely indicative of the activation barrier. The analytical fit has asymptotes at r_H–H = 0.770 Å and r_O–H = 0.903 Å.

Fig. 9 (a) Correlation of the distances r_O–H and r_H–H at the transition state, and the Heyrovsky activation energy Ω^‡_h for all dopants at various values of potential and hydrogen coverage. (b) Correlation of the promotion energy E_pd before hydrogen adsorption and Ω^‡_h at the corresponding potential and hydrogen coverage. The promotion energy captures the main trend of Ω^‡. Some particular features are exemplified by the Fe-system (horizontal), Pt-system (vertical), and Ti-system (gap crossing).

In terms of electronic structure, the previously discussed promotion energy E_pd captures the overall picture well. E_pd is positively correlated with adsorption energy, and relates to the Heyrovsky barrier as seen in Fig. 9b. Interestingly, the relation within each individual system does not necessarily follow the overall trend. This illustrates the different contributions to the activation mechanism (chemical adsorption energy as described by E_pd, and the electrode potential), and is especially visible for the Mn and Fe systems. These systems are hexagonally symmetric at PZC, but the ground state symmetry is spontaneously broken once a certain amount of excess negative charge is introduced. Additional charge further increases the distance between the dopant atom and the activated S, but during this transition the Fermi level (and hence potential) remains nearly constant, as the energies of unoccupied S-p states are lowered simultaneously with introduction of more electrons. In Fig. 9b this manifests as a local horizontal trend. For Fe and Mn, the changes in r_M–S are roughly 0.20 Å and 0.25 Å, respectively. Comparing with vertical trend systems such as Cu and Pt_θ=1/3, we find changes of roughly 0.015 Å. Within each system, E_pd is thus closely related to the distance r_M–S and does not seem to otherwise depend strongly on the amount of excess charge, unless it leads to the crossing of small gaps. Large change in E_pd due to gap crossing without significant geometrical change can be seen in Fig. 9b for the Sc and Ti dopants. In systems with small change in the promotion energy, the potential is the main contribution to barrier decrease.

4 Discussion

In all cases, single-atom doping results in a lowering of the PZC hydrogen adsorption energy. For the later 3d-metals, the doping results in breaking the local hexagonal symmetry, which leads to relatively larger affinity of the sulfur atom to hydrogen adsorption. This further leads to a wide range of HER-behavior for the dopants, from Volmer-limited to Heyrovsky-limited and even Tafel-limited in the case of Zn at moderate applied potential. Within the setting of the basal plane, particularly Co and Ni stand out with a good balance of the two barriers at moderate applied potential, though according to the kinetic model Ni is clearly more active at −0.5 V. At that point, the Ni-doped system exhibits barriers of 0.47 and 0.76 eV. In comparison, the Volmer and Tafel barriers on Pt(111) have with similar methodology been found to be 0.66 and 0.55 eV, respectively, at 0 V and a full monolayer H-coverage.⁴⁵ Therefore, the performance of the doped basal plane will be (at best) comparable to monocrystalline Pt(111) at a 0.5 V more negative applied potential. However, experiments often display comparable exchange currents at significantly lower potential offsets, see e.g. ref. 21. This coincides with the consensus that the pristine basal plane is not the main origin of MoS₂ activity, and edges and defects must be considered. When sites on certain facets are more active than others, the experimentally observed activity will depend largely on the morphology of prepared samples. Also the orientation of MoS₂ crystals with respect to the electrode substrate is important, as there is considerable resistance associated with electron transport between layers.⁵²

In the work by Humphrey et al.,²⁸ planar support is thought to have produced relatively low edge-content MoS₂ for the pristine case and with low levels of Co-doping. As noted there, the larger overpotential compared to other studies suggests that this activity is more representative of the basal plane. The distinction between the edges and basal-plane is illustrated directly in ref. 15, which reports respective overpotentials of around 0.3–0.4 V and 0.9–1.0 V. This inherent activity of the basal plane, while significantly lower than the edges, cannot be explained without introducing defects, as the pristine adsorption sites are much too unfavorable for evolution to occur. Most synthesized MoS₂ does however contain a significant number of sulfur vacancies, with certain deposition methods yielding stoichiometries of MoS_1.6⁵³ and MoS_1.8,²⁵ although this total deficiency will also depend in part on the stoichiometry and prevalence of the edge terminations. These vacancies create local sites with adsorption energies of roughly 0.1 eV. However, neighboring sites are not significantly activated and (assuming evenly dispersed S-vacancies) the reaction would still be limited to a Volmer–Heyrovsky mechanism. More complex defect configurations could possibly facilitate the Tafel mechanism. For reference, the Tafel slope of the basal plane of pristine MoS₂ under acidic conditions has been measured to be around 120 mV dec⁻¹,⁵⁴ indicative (but not conclusive) of a mechanism rate-determined by the Volmer step, from which neither the Tafel nor Heyrovsky path can be disregarded.

Importantly, the aforementioned study finds that low levels of Co-doping is detrimental to the inherent activity, in contradiction to the results of the present study in which only stoichiometric MoS₂ is considered. This indicates that the combination of atomic doping and intrinsic defects can lead to overall deactivation of the basal plane. In ref. 28, this is supported by the calculated adsorption energies and equilibrium defect configurations, and it would be interesting to study this problem in terms of activation barriers and possible reaction paths of the Co- and Ni-doped defect systems.

We summarize that a direct comparison of our results with experiment is not possible without also considering MoS₂ edge sites and defects in the basal plane, as well as other possible doping configurations and interplay with inherent sulfur vacancies. A systematic investigation of activation barriers and coverage dependence across these systems is needed for drawing solid conclusions.

The overall trend of relative activity is clear in Fig. 8b, the optimum lies in the vicinity of ΔΩ_H ≈ 0 eV, though we note that the peaks of the directly fitted volcanoes are in this case centered at 0.3 eV. The discrepancy between the solid line and the points illustrates the expected inaccuracies associated with ΔG_H as a descriptor, as the relation to the activation energy is only implicit. Overall the description works well, but comparing individual systems in terms of only ΔG_H will be prone to errors. On a larger scale, it is important to note that the condition ΔG_H ≈ 0 eV is not sufficient to guarantee high absolute HER activity, but is merely a necessary condition to be near the optimum for a given class of systems. As an example, we see in this study that the Heyrovsky mechanism on local impurity-induced S-sites in MoS₂ behaves markedly differently from e.g. the Tafel mechanism on uniformly covered Pt(111), even if the H adsorption energy is near zero in both cases.

5 Conclusions

Methodology for computing reaction and activation energies under electrochemical conditions corresponding to a specified applied voltage was used to investigate the HER mechanism on the basal plane of 2H-MoS₂ with single-atom 3d-metal doping. The effect of the aqueous solvent is included by using a cluster of a four H₂O molecules and a proton (Eigen cation) in the neighborhood of the active site and then a polarizable continuum solvent for the rest of the solvent phase. The effect of the various transition–metal dopants spans a wide range of activation energies and within this scope Ni stands out with the best overall activity at moderate negative applied potential. The reaction barriers can be correlated with the potential-dependent promotion energy E_pd, and the modelled kinetics has the characteristic volcano-shape with respect to the hydrogen adsorption energy. The results were compared with the commonly used ΔG_H model, and the condition ΔG_H ≈ 0 eV was shown not to be sufficient for predicting a low reaction barrier in these systems.

The calculations used a single monolayer of 2H-MoS₂, but the results are expected to be representative for multi-layer slabs as the adsorption energies are found to be similar. The reaction was found to proceed predominantly through a Volmer–Heyrovsky path, where the dopant-activated sulfur sites provided a large reduction in the adsorption energy, and therefore, also in the activation energy of the Volmer step which on pristine MoS₂ is so unfavorable that it essentially does not occur. Thus, the inherent activity of MoS₂ observed in experiments must come from defect sites (e.g. sulfur vacancies on the basal plane) and/or from edge sites. These sites are expected to be more active than the pristine basal plane, and a full picture of the various possible sites is necessary for a more in-depth comparison with experimental data. Such an investigation is planned for future work.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We thank K. Laasonen for discussions. The calculations were performed on resources provided by Sigma2 – the National Infrastructure for High Performance Computing and Data Storage in Norway, project No. NN9497K. J. A. acknowledges financial support from the Academy of Finland, project No. 322832 “NANOIONICS”. H. J. acknowledges financial support from the Icelandic Research Fund, project No. 207283-053.

Notes and references

1. F. Dawood, M. Anda and G. Shafiullah, Int. J. Hydrogen Energy, 2020, 45, 3847–3869.
2. J. Holladay, J. Hu, D. King and Y. Wang, Catal. Today, 2009, 139, 244–260.
3. M. S. Faber and S. Jin, Energy Environ. Sci., 2014, 7, 3519–3542.
4. Y. Pan, C. Zhang, Y. Lin, Z. Liu, M. Wang and C. Chen, Sci. China: Mater., 2020, 63, 921–948.
5. H. Jin, C. Guo, X. Liu, J. Liu, A. Vasileff, Y. Jiao, Y. Zheng and S.-Z. Qiao, Chem. Rev., 2018, 118, 6337–6408.
6. C. G. Morales-Guio, L.-A. Stern and X. Hu, Chem. Soc. Rev., 2014, 43, 6555–6569.
7. J. Xie, X. Yang and Y. Xie, Nanoscale, 2020, 12, 4283–4294.
8. J. Xie, J. Qi, F. Lei and Y. Xie, Chem. Commun., 2020, 56, 11910–11930.
9. B. Hinnemann, P. G. Moses, J. Bonde, K. P. Jørgensen, J. H. Nielsen, S. Horch, I. Chorkendorff and J. K. Nørskov, J. Am. Chem. Soc., 2005, 127, 5308–5309.
10. J. Bonde, P. G. Moses, T. F. Jaramillo, J. K. Nørskov and I. Chorkendorff, Faraday Discuss., 2009, 140, 219–231.
11. J. D. Benck, Z. Chen, L. Y. Kuritzky, A. J. Forman and T. F. Jaramillo, ACS Catal., 2012, 2, 1916–1923.
12. J. Kibsgaard, Z. Chen, B. N. Reinecke and T. F. Jaramillo, Nat. Mater., 2012, 11, 963–969.
13. J. Xie, J. Zhang, S. Li, F. Grote, X. Zhang, H. Zhang, R. Wang, Y. Lei, B. Pan and Y. Xie, J. Am. Chem. Soc., 2013, 135, 17881–17888.
14. T. F. Jaramillo, K. P. Jørgensen, J. Bonde, J. H. Nielsen, S. Horch and I. Chorkendorff, Science, 2007, 317, 100–102.
15. S. M. Tan, A. Ambrosi, Z. Sofer, Š. Huber, D. Sedmidubský and M. Pumera, Chem. – Eur. J., 2015, 21, 7170–7178.
16. D. Kong, H. Wang, J. J. Cha, M. Pasta, K. J. Koski, J. Yao and Y. Cui, Nano Lett., 2013, 13, 1341–1347.
17. J. Xie, H. Zhang, S. Li, R. Wang, X. Sun, M. Zhou, J. Zhou, X. W. D. Lou and Y. Xie, Adv. Mater., 2013, 25, 5807–5813.
18. J. Xie, H. Qu, J. Xin, X. Zhang, G. Cui, X. Zhang, J. Bao, B. Tang and Y. Xie, Nano Res., 2017, 10, 1178–1188.
19. H. Li, C. Tsai, A. L. Koh, L. Cai, A. W. Contryman, A. H. Fragapane, J. Zhao, H. S. Han, H. C. Manoharan, F. Abild-Pedersen, J. K. Nørskov and X. Zheng, Nat. Mater., 2016, 15, 48–53.
20. T. Zhang, H. Zhu, C. Guo, S. Cao, C.-M. L. Wu, Z. Wang and X. Lu, Catal. Sci. Technol., 2020, 10, 458–465.
21. J. Deng, H. Li, J. Xiao, Y. Tu, D. Deng, H. Yang, H. Tian, J. Li, P. Ren and X. Bao, Energy Environ. Sci., 2015, 8, 1594–1601.
22. K. Sun, L. Zeng, S. Liu, L. Zhao, H. Zhu, J. Zhao, Z. Liu, D. Cao, Y. Hou, Y. Liu, Y. Pan and C. Liu, Nano Energy, 2019, 58, 862–869.
23. H. Wang, C. Tsai, D. Kong, K. Chan, F. Abild-Pedersen, J. K. Nørskov and Y. Cui, Nano Res., 2015, 8, 566–575.
24. D. Merki, H. Vrubel, L. Rovelli, S. Fierro and X. Hu, Chem. Sci., 2012, 3, 2515–2525.
25. D. Escalera-López, Y. Niu, J. Yin, K. Cooke, N. V. Rees and R. E. Palmer, ACS Catal., 2016, 6, 6008–6017.
26. Y. Shi, Y. Zhou, D.-R. Yang, W.-X. Xu, C. Wang, F.-B. Wang, J.-J. Xu, X.-H. Xia and H.-Y. Chen, J. Am. Chem. Soc., 2017, 139, 15479–15485.
27. T. H. M. Lau, X. Lu, J. Kulhavý, S. Wu, L. Lu, T.-S. Wu, R. Kato, J. S. Foord, Y.-L. Soo, K. Suenaga and S. C. E. Tsang, Chem. Sci., 2018, 9, 4769–4776.
28. J. J. L. Humphrey, R. Kronberg, R. Cai, K. Laasonen, R. E. Palmer and A. J. Wain, Nanoscale, 2020, 12, 4459–4472.
29. M. Hakala, R. Kronberg and K. Laasonen, Sci. Rep., 2017, 7, 15243.
30. B. Hammer, L. B. Hansen and J. K. Nørskov, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 7413–7421.
31. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
32. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
33. B. Himmetoglu, A. Floris, S. de Gironcoli and M. Cococcioni, Int. J. Quantum Chem., 2014, 114, 14–49.
34. J. P. Perdew, M. Ernzerhof and K. Burke, J. Chem. Phys., 1996, 105, 9982–9985.
35. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158–6170.
36. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
37. G. Henkelman, B. P. Uberuaga and H. Jónsson, J. Chem. Phys., 2000, 113, 9901–9904.
38. G. Henkelman and H. Jónsson, J. Chem. Phys., 1999, 111, 7010–7022.
39. S. Trasatti, J. Electroanal. Chem. Interfacial Electrochem., 1972, 39, 163–184.
40. J. K. Nørskov, T. Bligaard, A. Logadottir, J. R. Kitchin, J. G. Chen, S. Pandelov and U. Stimming, J. Electrochem. Soc., 2005, 152, J23.
41. J. K. Nørskov, J. Rossmeisl, A. Logadottir, L. Lindqvist, J. R. Kitchin, T. Bligaard and H. Jónsson, J. Phys. Chem. B, 2004, 108, 17886–17892.
42. K. Mathew, R. Sundararaman, K. Letchworth-Weaver, T. A. Arias and R. G. Hennig, J. Chem. Phys., 2014, 140, 084106.
43. K. Mathew, V. S. C. Kolluru, S. Mula, S. N. Steinmann and R. G. Hennig, J. Chem. Phys., 2019, 151, 234101.
44. J. Rossmeisl, E. Skúlason, M. E. Björketun, V. Tripkovic and J. K. Nørskov, Chem. Phys. Lett., 2008, 466, 68–71.
45. M. Van den Bossche, E. Skúlason, C. Rose-Petruck and H. Jónsson, J. Phys. Chem. C, 2019, 123, 4116–4124.
46. S. Trasatti, J. Electroanal. Chem. Int. Electrochem., 1986, 209, 417–428.
47. E. Santos and W. Schmickler, Chem. Phys. Lett., 2004, 400, 26–29.
48. N. Wakabayashi, H. G. Smith and R. M. Nicklow, Phys. Rev. B: Condens. Matter Mater. Phys., 1975, 12, 659–663.
49. T. Böker, R. Severin, A. Müller, C. Janowitz, R. Manzke, D. Voß, P. Krüger, A. Mazur and J. Pollmann, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 64, 235305.
50. K. F. Mak, C. Lee, J. Hone, J. Shan and T. F. Heinz, Phys. Rev. Lett., 2010, 105, 136805.
51. M. Liu, M. S. Hybertsen and Q. Wu, Angew. Chem., Int. Ed., 2020, 59, 14835–14841.
52. Y. Yu, S.-Y. Huang, Y. Li, S. N. Steinmann, W. Yang and L. Cao, Nano Lett., 2014, 14, 553–558.
53. M. J. Cuddy, K. P. Arkill, Z. W. Wang, H.-P. Komsa, A. V. Krasheninnikov and R. E. Palmer, Nanoscale, 2014, 6, 12463–12469.
54. C. L. Bentley, M. Kang, F. M. Maddar, F. Li, M. Walker, J. Zhang and P. R. Unwin, Chem. Sci., 2017, 8, 6583–6593.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3cp00516j

---