How to Dope the Basal Plane of 2H-MoS 2 to Boost the Hydrogen Evolution Reaction?
Molybdenum disulfide (MoS 2 ) is considered one of the most likely materials that could be turned into low-cost hydrogen evolution reaction (HER) catalysts to replace noble metals in acidic medium. However, several challenges prevent MoS 2 from being truly applicable, including limited number of active sites (typically only the edges are active) and poor conductivity. In this work, we perform an extensive density functional theory (DFT) screening of substitutional doping as a possibility to activate the otherwise inert basal surface. We assess 17 Earth abundant elements for molybdenum doping and 5 elements (N, O, P, Se and Te) for sulfur substitution. Systematically determining the preference of the metallic dopants to be located on the edges rather than in the basal plane, we reveal that most dopants are much more likely to be incorporated at the edges, suggesting that advanced synthesis methods are required to obtain basal-plane doped catalysts. The latter may, however, feature many more active sites per MoS 2 formula unit, motivating our study on the properties of such substitutionally doped surfaces. For the first time for such a screening study, we explore not only the 1 adsorption of H, but also of OH and H 2 O to explore the solvent effect since the reaction takes place in an aqueous medium. Two additional phenomena that could hinder the hydrogen production at these sites are investigated, namely H 2 S release and the (local) segregation/dispersion tendency of the dopants in the basal surface. Moreover, to assess the electrocatalytic activity, we take the electrochemical potential explicitly into account via grand canonical DFT in combinations. Compared with pristine MoS 2 nanosheets, our results show that most doping elements significantly enhanced the electrocatalytic activity. Considering all assessed factors, we identify the most promising systems: Dimers of Ti, Zr and Hf and the substitution of S by P are predicted to lead to stable active sites on the basal plane with overpotentials of about 0.2 V.
Introduction<!>Computational details<!>Results and Discussion<!>Relative Stability of<!>MoS 2 Basal Plane Doping<!>Monomer:<!>Conclusion<!>Acknowledgement
<p>Hydrogen production by water electrolysis is an intensely studied process for the mediumterm storage of intermittent, decarbonized electricity such as solar, wind and tidal energy. 1,2 Furthermore, "clean" hydrogen could also become a widely used energy vector for heavy transportation, replacing fossil fuels. However, the high cost of platinum electrodes represents one of the barriers towards the large-scale implementation of electrolytic generation of H 2 with the Proton Exchange Membrane (PEM) technology. 3,4 In recent years, transition metal dichalcogenides (TMDs) with their intriguing structural and electronic properties have attracted increasing interest. [5][6][7] Their large-scale synthesis, based on the processes of liquid exfoliation, chemical vapor deposition and especially impregnation followed by sulfidation, 8 facilitates their use in various fields. 9,10 The most important TMD is MoS 2 , which has attracted great interest as a low-cost and environmentally friendly electrocatalyst for HER mechanism. [11][12][13] Furthermore, MoS 2 catalysts are already produced on large scales for hydro-treatment and hydro-cracking reactions for oil refining.</p><p>In general, the activity of the catalyst depends mainly on the total number of active sites and on their intrinsic activity. For pristine MoS 2 , both theoretical calculations and subsequent experimental results have shown that promising sites are generally confined to the small region of exposed edges. 14 The predominantly exposed basal plane of the semiconducting 2H-MoS 2 phase is, however, catalytically inactive in the absence of defects. 13,15 Thus, the search for a method to activate the basal plane to exploit its large accessible surface area had led to extensive research. One of the many attempts focuses on the phase transformation of 2H to 1T polytype, which has promising sites on its basal plane but is metastable. [16][17][18] It was also reported that the 2H basal plane could be activated by the introduction of various types of defects and in particular sulfur vacancies. 19,20 The conductivity of MoS 2 can also be improved by adding conductive supports such as metals and carbon materials. 21,22 A common strategy to activate MoS 2 is doping. Doping is a typical strategy to precisely control the electronic, 23 optical, 24 magnetic, 25 and catalytic properties of 2D materials. 26 Four different doping strategies with different effects have to be distinguished. First, surface charge doping has attracted attention for modulating the electronic behavior of semiconductors. 27 It takes advantage of the interaction between the host material and molecules, ions and solids. If the dopant is in contact with the host material a charge-transfer takes place: oxidizing molecules are p-type dopants, 28,29 while reducing molecules are n-type dopants. [30][31][32] Second, small atoms can lead to intercalation compounds. 33 For MoS 2 , the intercalation of atoms that dissociate into small cations such as H + plus an electron leads to lattice distortion and a significant change in electronic structure due to electron donation. This approach allows to stabilize the 1T-like polymorph. 34 In addition to light elements , heavy elements such as Pt, Ru, Au and Pd have also been successfully intercalated in MoS 2 by Chen et al. conductivity and high number of active sites. Very recently, based on DFT computations, Ti doping was also suggested to be promising for increasing the HER activity. 49 As an alternative, also the substitutional doping of sulfur, rather than Mo, has been tested. For example, Ye et al. showed the improvement of the catalytic activity of MoS 2 by preparing MoS 2 − xP x solid solutions. 50 Similarly, N was theoretically and experimentally found to tune the HER activity of MoS 2 both by generating additional active sites and by increasing the conductivity. 51,52 It is to be noted that phosphorus doping might also "accidentally" originate from phosphate additives during the synthesis process of MoS 2 . 53 Isoelectronic dopants (such as Se and W) have also been tested. They are easier to incorporate but do not significantly change the electronic structure, so that only morphology related properties of MoS 2 are tuned, typically by increasing the number of active edge or basal plane defect sites. 54,55 All these reports are stimulating further engineering MoS 2 via dopants to achieve high HER performance. However, synthetic issues (e.g., non-scalable synthesis, unwanted phases...) and lacking long-term stability have, so far, prevented the development of an MoS 2 -based HER catalyst to be applicable on large scale 56 Therefore, identifying and understanding the impact of various dopants on the HER activity on MoS 2 remains an active topic.</p><p>In this work, we explore the influence on the HER acitivty of the basal plane of 2H-MoS 2 of substitutional metal-doping (17 non-noble metals) and/or non-metal doping (N, O, P, Se, Te) via an in silico screening. The primary criterion to classify them is the hydrogen adsorption energy ∆G H , which should be close to 0 eV. 57 If ∆G H is very negative, H will be strongly adsorbed and poison the catalyst surface, and if it is very positive, the Volmer step (adsorption of H + + e -) might never occur. We consistently assess the tendency of the dopants to cluster together or to form single sites and determine the influence of this tendency on the catalytic properties. After an initial screening, we explicitly consider the electrochemical potential to study the reaction mechanism on promising sites, including the reactivity with water. This assessment, together with the determination of the stability of the doped basal planes with respect to H 2 S-loss, and the impact of the dopants on the conductivity, provide a comprehensive overview on the prospects to induce catalytic activity on the defect free MoS 2 basal plane via substitutional doping.</p><!><p>Density functional theory (DFT) computations were performed within the Vienna ab initio simulation package (VASP). 58 Electron-ion interactions were modeled by projector augmented wave potentials. 59,60 The generalized gradient approximation in the flavor of the dispersion corrected Perdew-Burke-Ernzerhof functional (PBE-dDsC) was applied to approximate exchange and correlation effects. [61][62][63] For the plane-wave basis set expansion, the cut-off energy was set to 500 eV. The precision setting of VASP is set to "accurate". The Brillouin zone is integrated by a Monkhorst k-points mesh. For the description of the (doped) p(4 × 4) cell of the (0001) facet only the Γ-point was considered. For the edges, we adopted a 3 × 3 × 1 k-point mesh. The convergence criterion for the self-consistency process is set to 10 −5 eV for the optimization of the wave function. The maximum forces are converged to 0.025 eV during geometry optimization. Approximately 15 Å of vacuum is added both above and below the slabs in the z-direction. 63 For studying substitutional doping, we have assessed the energy difference between locating the dopant on the Mo-or S-edge compared to placing it in the basal plane. This is carried out with a fully relaxed symmetric six-layer p(1×2) edge model (see Figure 1).</p><p>Our first stage in the screening of substitutional doping of the basal plane consisted either the replacement of one sulfur atom by Y= N, O, P, Se, and Te or by replacing a triangle of molybdenum atoms by X= Ti, V, Cr, Mn, Fe, Co, Ni, Zn, Ga, Zr, Nb, Sn, Sb, Hf, Ta, W and Bi. For X=[Cr, Mn, Fe, Co, Ni] we have considered spin polarization and tested various initial spins to identify the most stable magnetic state. For certain X and Y, we have also investigated a mixed doping denoted by X/Y (see Figure 2 for a schematic summary). For  this first stage of the screening, the model system consists of 3 MoS 2 layers, with the lowest layer fixed in its bulk position and the dopants located in the top layer. We tried three doping densities (θ) defined as follows:</p><p>Given the system size (Mo 48 S 96 ), the doping density corresponds to to 0.69 %, 1.39 % and 2 % for 1X/S, 2X/S, and 3X/S, while only 0.69 % is tested for Mo/Y.</p><p>An important parameter is the tendency of dopants to either be dispersed in MoS 2 or to cluster together. To quantify this tendency, we determine the energy to break up dimers and trimers, ∆G D and ∆G T , respectively. Taking the fully dispersed model as a reference, we define:</p><p>and</p><p>where X is the doping element.</p><p>The catalytic activity of all the mentioned structures is assessed by the hydrogen adsorption energy. For the basal planes, hydrogen is adsorbed on the "bent" position in agreement with previous studies. 20 The hydrogen adsorption energy (∆G H ) is computed using the following equation:</p><p>where E H and E slab are the total energies of the system with and without hydrogen adsorbed on the active site, respectively, and G H 2 is the total energy of the hydrogen molecule:</p><p>-TS where T is the room temperature (298.15 K) and S is the rotational and translational entropy at standard pressure.</p><p>Considering ∆G H and the relative stability of the dimers and trimers, we selected the promising structures to perform grand-canonical DFT computations. For these computations we rely on symmetric five layer structures (Mo 80 S 160 ), where the central layer is frozen in its bulk geometry. For these systems, the doping densities are equal to 0.83 %, 1.67 % and 2.5</p><p>% for Mo and 0.83 % for S doping. These symmetric systems have been utilized to determine the potential dependent adsorption energy of H, OH and H 2 O. For these grand-canonical DFT computations, the linearized Poisson-Boltzmann model as implemented in VASPsol is applied. 64 In VASPsol, the interface is defined as a function of the electron density, which is compared to a "critical" density, ρ c . In line with our previous work 20 we set the iso-density value to to 0.00025 e − / Å−1 . This is the same value as previously also used in the context of solvated alkali-metal ions. 65,66 The Debey screening length is set to 3 Å, corresponding to a 1 M electrolyte. The potential dependence is studied by adapting the surface charge.</p><p>Practically, the charge of the system is varied from −1 to +1 in steps of 0.2 e − . The acquired connection between the grand-canonical energy and the electrochemical potential is then fitted to an analytical formula, which is a parabola for conductors and adapted to semi-conductors by accounting for the band-gap. 20,67,68 For the grand-canonical DFT results, the adsorption energy ∆G ads of the different species are calculated according to:</p><p>where E tot (U ) and E slab (U ) are the potential-dependent energies of the symmetric MoS 2 surface with and without adsorbates, respectively. If the adsorbate is hydrogen (i.e., one hydrogen atom per surface), the reactant (R) is defined as:</p><p>where G H 2 (g) is the free energy of the hydrogen molecule, e is the elementary charge, and U is the applied potential with respect to the standard hydrogen electrode. This approach relies on the seminal paper introducing the computational hydrogen electrode, where Norskov and coworkers pointed out that the chemical potential of protons and electrons can be conveniently replaced by the one of a hydrogen molecule and a simple energy correction to account for the electrode potential. 69 Eq. 6 requires, in general, a pH dependence. Since we here study HER in acidic conditions, we assume a pH of 0 for which Eq. 6 is valid.</p><p>For the adsorption of the water molecule, R = 2G</p><p>, so that it formally should be considered as an adsorbed hydroxyl "radical" and not to be confused with adsorbed hydroxide "ion" OH -. The adsorbed *OH can be an intermediate in HER (see Fig. 3), as long as its formation is exothermic and the corresponding reaction energies</p><p>are close to zero.</p><p>We also compute the electrical conductivity (σ) of the most promising doped systems using BoltzTraP2, 71 which solves the semiclassical Boltzmann transport equation under the constant relaxation time approximation. We obtain the electrical conductivity at 300 K around the Fermi energy level. This allows us to obtain a numerical estimate of the changes in conductivity upon doping.</p><!><p>The adsorption of H on the pristine basal plane of MoS 2 is unstable with a positive ∆G H of up to 1.67 eV, while an ideal ∆G H should be close to zero. We have explored two strategies to activate the inert basal plane: Curvature of MoS 2 , modelled by nanotubes, and substitutional doping. However, since the activity of MoS 2 nanotubes is predicted to be disappointing, these results are only discussed in the SI (see section S1). Figure S3 represents the detailed results for all dopants and positions tested, while Figure 4 provides a summary of the main findings. All dopants have an energetic preference for being located at the edges. However, the extent of this preference strongly depends on the particular metal. Most elements are most stably located on the 50%S S-edge and only Nb, Zr, Hf and Ta are better stabilized on the 50%S Mo-edge. Nevertheless, experimentally it is well documented that many elements can be doped into the basal plane thanks to experimental techniques such as one-step sintering 52 and single-step hydrothermal reactions 72 leading to a good control of the amount and the position of the doped atoms. Hence, we have chosen a "critical" energy difference of 2 eV with respect to the most stable configuration: below this energy we expect that the substitutional doping should be feasible, while above the basal-plane doping seems largely compromised. We observe that only V, Nb, Ta, Cr and W are able to dope the MoS 2 basal plane with ∆E lower than 2 eV. Ti, Zr and Hf are also located close to this threshold and could, therefore, probably be doped in the basal plane as well. All other elements are strongly stabilized at the edges compared to the basal plane.</p><!><p>In view of the experimental reports of basal-plane doped MoS 2 48,50,51,73 and the observation that the intrinsic activity of the edges is hard to improve 47 (i.e., edge-engineering mainly increases the number of active sites, not their thermodynamic overpotential), we investigate the large pannel of non-noble metal dopants for activating the basal plane in the next section.</p><!><p>The geometry optimization of doped MoS 2 leads in most cases to a visible deformation. The most common deformation observed is a slight deviation from planarity as exemplified in Figure 5 a. In some cases, especially for large dopants (Bi, P, Sb) one of the sulfur atoms finds itself practically expelled from the surface (Fig. 5 b, c), while in the case of 3Ni/S a S-S dimer is formed (Fig. 5 d). Note that these geometries and others with adsorbed hydrogen are publicly available and additional geometries shown in the SI (see Fig. S5).</p><p>To select the most promising dopants, Computational Hydrogen Electrode (CHE)-based DFT calculations of the hydrogen adsorption energies are performed. Only for the promising dopants we perform the grand-canonical DFT computations for H, OH and H 2 O adsorption.</p><p>The influence of the number of substituted Mo was investigated by a systematic series of systems consisting in single-atom doping, dimer doping and a triangular trimer doping. The CHE results led to the exclusion of Cr and W, since ∆G H was found to be very positive at all doping densities tested (the lowest ∆G H value for Cr and W is 1.11 eV and 1.91 eV, respectively). For the other elements, the HER activity is significantly improved compared to the pristine basal plane. We distinguish three families depending on their interaction with H. The first features ∆G H closest to 0 eV when θ= 2% (3Mn/S, 3Ni/S, 3V/S, 3Nb/S, 3Sn/S and 3Ta/S), the second when θ=1.39% (2Ti/S, 2Zr/S, 2Hf/S, 2Sb/S and 2Bi/S) and the third is when θ= 0.69% (1Co/S, 1Fe/S, 1Zn/S and 1Ga/S) (see Table S3). We also add to the list 3Ti/S, 3Zr/S and 3Hf/S as ∆G H is very close to 0 eV for these systems (-0.12 eV, 0.17 eV and 0.2 eV respectively). Since the hydrogen adsorption energy strongly depends on the doping density, we determine the energy required to break up dimers and trimers (see Eq. 2 and 3, respectively) for all fourteen relevant dopants. Figure 6 depicts these reaction energies, with the single-atom dopant as the reference (∆G M = 0, M: monoatomic). Lines connecting the energies have no physical meaning, but are added as a guide to the eye to facilitate the interpretation, i.e., whether or not a given dopant will tend to form "clusters" or is more stable as dispersed single-atoms in the basal plane.</p><!><p>Dimer:</p><p>Trimer: When ∆G is negative, the dopant is found as a single atom in the basal plane, while a positive reaction indicates that cluster formation is thermodynamically likely. Three dopants show a uniform negative slope, indicating that single-atom doped MoS 2 is thermodynamically favored for V, Nb and Ta. Titanium, hafnium and zirconium show a weak preference for dimer formation. All the other dopants show a more (Ni) or less (Mn) pronounced tendency to cluster into trimers. This segregation tendency could indicate that metal-segregated hybrid structures would be obtained experimentally with most synthesis techniques. Thus, some of the structures mentioned above to be promising for HER according to ∆G H are not favorable to be formed based on the results showed in Figure 6, in particular 3V/S, 3Nb/S, 3Ta/S, 2Sb/S, 2Bi/S, 1Co/S, 1Fe/S, 1Zn/S and 1Ga/S. In view of this additional "constraint" we have revised the list of systems that warrant a more detailed (grand-canonical DFT) investigation. In particular, we have added 3Zn/S, 3Bi/S, 2Sb/S, 2Ga/S, 1Nb/S and 1V/S.</p><p>The second type of doping is to substitute the sulfur atom. The tested elements are N, O, P, Se, and Te (see Table S4). For O, Se and Te, the adsorption energies of hydrogen are 1.1 eV, 2.0 eV and 1.9 eV, respectively. Hence, they clearly do not activate the basal plane.</p><p>Nitrogen and phosphorous, however, gave promising results with ∆G H below 0 eV (-1.5 eV and -0.7 eV, respectively). As ∆G H is very negative, we determined the adsorption energy for a second hydrogen atom. However, this second hydrogen atom is not bound (adsorption energies of 2.5 eV and 1.5 eV for Mo/N and Mo/P, respectively), in line with expected valencies of these group 15 elements. Despite the strong H adsorption, we will assess the interaction of these sites with water to determine if this might lead to HER active sites in analogy to some anti-site defects studied before. 20 Finally, the mixed substitutional doping has been extensively explored. However, only few combinations turned out to be more interesting than the (conceptually simpler) single doped systems (see Table S5). In particular, 3Bi/P, 3Nb/P, 2Nb/P and 1Nb/P have been identified as sufficiently promising to be included in the more detailed grand-canonical DFT study below.</p><p>The possible mechanisms to produce H 2 at the promising active sites are Volmer-Heyrovsky if ∆G H is close to 0 eV or Volmer-Tafel if ∆G H is negative and a second hydrogen can be adsorbed. This second scenario is encountered for (3Zn, 3Sb, 3Bi, 3Ga, and 2Bi). How-ever, upon the (exothermic) adsorption of the second hydrogen atom, all these cases lead to the release of H 2 S. As a consequence an S-defect is created. While such defects might be catalytically active, they are beyond the scope of the present investigation that focuses on activating the basal plane in the absence of defects. Given that HER is performed in aqueous electrolytes, we also assess the interaction of the potential active sites with H 2 O an electrochemically formed OH. These species can either only weakly interact and thus be of no consequence for HER or they can "block" HER, typically forming stable surface OH* species or participate in HER, if oxygen is part of the active site and O-H bonds are easily broken and formed.</p><p>The influence of the electrochemical potential on the HER activity of the selected systems Having identified the promising structures, we proceed to evaluate their HER activity by considering the influence of the electrochemical potential and the presence of other species, OH* and H 2 O, since the presence of water at the electrochemical interface can have profound effects on the activity beyond the macroscopic effects captured by the Poisson-Boltzmann equation. 20 In the supporting information (see section S3) we display and discuss typical adsorption energy behaviors of H, OH and H 2 O as a function of the electrochemical potential.</p><p>Figure 7 reports the relevant adsorption energies at the thermodynamical equilibrium potential for HER, i.e., 0 V vs SHE. Note, that these energies differ from the values for neutral surfaces used for the prescreening, since we are applying a grand-canonical DFT framework. The explicit inclusion of the electrochemical potential leads to an adjustment of the CHE results discussed above. Thus, we arrive at a new classification based on this more reliable level of theory. First, several systems (1V/S, 3Mn/S, 3Nb/S, 1Nb/S, 1Bi/S and 3Ta/S) feature ∆G H (U ) equal to or greater than 0.5 eV and are, thus no longer relevant for the typical HER mechanism (Fig. 7a). Since these sites seem unreactive with repsect to H*, we did not consider OH* which tends to be less strongly adsorbed. Second, for 3Zn/S, 3Sb/S, 3Bi/S and 3Ga/S ∆G H (U) is very negative (≤ −0.5 eV; see Fig 7b). However, these systems still suffer from the instability with respect to H 2 S release. The last group of subsitutional MO doping is when ∆G H (U) is close to 0 eV (Fig. 7c). This promising class of doped systems includes 3Ti/S, 2Ti/S, 3Ni/S, 1Co/S, 1Fe/S, 3V/S, 1Zn/S, 3Zr/S, 2Zr/S, 3Sn/S, 3Hf/S, 2Hf/S, 2Sb/S, 2Bi/S, and 1Ga/S. For these active sites the Volmer step is likely to occur at moderate over-potentials and, since the bonding is not strong, H 2 can be subsequently produced via a Heyrovsky step. However, these active sites might be blocked by OH*. As shown in Figure 7, we find that ∆G OH (U) is positive for these active sites. Hence, OH* is not expected to block them. Furthermore, H 2 O only interacts weakly with these surfaces, confirming the absence of strong interference of the solvent with HER, except for 3Ga/S where ∆G H 2 O (U) =-0.55 eV. In this particular case, the deformation of the surface exposes on of the Ga dopants, so that it can interact directly with the oxygen of the water molecule. We conclude that for these systems the promising sites are not blocked and the Volmer-Heyrovsky mechanism is possible at moderate overpotentials.</p><p>We now turn to the second type of doping, in which the sulfur atom is substituted, including mixed doping. We adsorb H, OH and H 2 O. Figure 7d, shows that water adsorption remains thermoneutral. Hydrogen adsorption, however, is strongly exothermic for most cases. Niobium co-doped with phosphorus (1Nb/P, 2Nb/P and 3Nb/P) stabilizes a second adsorbed hydrogen, although no strongly enough to be relevant for HER. On the other hand, ∆G H (0 V) is very close to 0 eV for 3Bi/P. Thus, the Volmer-Heyrovsky mechanism could occur on this particular active site. Similar to hydrogen adsorption, ∆G OH (0 V) is negative for all systems except Mo/N. Consequently, OH* might block most of these sites, except when Volmer-Heyrovsky-type reaction mechanism involving OH* and H 2 O (see Figure 3) are favorable, which we explore in more detail in the next section.  To assess the occurrence of H 2 O/ OH cycles to generate H 2 , we calculate ∆G 1 and ∆G 2 according to Eq. 7 and Eq. 8 as a function of the electrochemical potential. Only for three systems out of all the doped structures the formation of OH* is feasible and ∆G 1 and ∆G 2 are close to or smaller than 0 eV at overpotentials ≤ 0.5 V, namely Mo/P, 3Bi/S, and 3Bi/P with thermodynamic overpotentials η T D equal to 0.2 V, 0.3 V and 0.5 V, respectively. Figure 8 illustrates the active site with the lowest thermodynamic overpotential in the series (0.2 V), namely Mo/P. Incidentally, this system is the one that is most likely to be achievable experimentally among the screened S-substituted systems. 50,74 Conductivity and identification of most promising systems Grand-canonical DFT reveals the electronic structure of the various systems as either effectively metallic (e.g., H@2Zr/S following a parabola fit) or as semi-conducting (e.g., H@Mo/P which shows a piece-wise parabola on the reductive and oxidative side of the fundamental gap) as shown in Figure S4. However, this indirect probe of the electronic structure remains very qualitative. In order to quantitatively estimate the conductivity of the promising doped systems, we resorted to the semi-classical transport equations as implemented in BoltzTraP2.</p><p>The obtained electrical conductivities are summarized in the last column of Figure S7. Most of these doped systems are conductive, the exceptions being Mo/P and 3Ni/S. Hence, according to our computations, doping of the basal plane not only leads to novel active sites, but also lifts one of the major limitations of MoS 2 electrocatalysts, namely their low electrical conductivity.</p><p>Figure S7 summarizes the information from the various tests performed, including relative energies of placing dopants in the basal plane or at the edges and the tendency (or not) of dopants to cluster together. The release of H 2 S and the geometrical distortion, especially upon H adsorption, are also summarized. Furthermore, all the relevant thermodynamic overpotentials, including at the OH* site, are reported. Based on these various criteria, we arrive at our final ranking: First of all, the considered site needs to produce H 2 at potentials ≥ -0.2 V vs SHE. Second, the corresponding mono-/dimer-/trimer doped system should be stable with respect to clustering/dispersion. Finally, doping of the basal plane compared to doping of the edge site should be relatively easy. From all the Mo substitutionally doped systems tested, only 2Ti/S, 2Zr/S and 2Hf/S are found to be truly promising, i.e., satisfying all three criteria. These sites generate H 2 via the Volmer-Heyrovsky mechanism. Mo/P is also considered promising, as it releases H 2 via the OH/H 2 O cycle and with a potential as low as 0.2 V. Given that ZrS 2 and HfS 2 have been reported as reasonable HER catalysts 75 and experiments indicate that doping with titanium and zirconium is possible in MoS 2 , 76,77 we expect that our predictions will spur corresponding experiments. Similarly, substitutional P doping of MoS 2 seems experimentally feasible 74 and our computations suggest that OH* species might not only be responsible for the experimentally observed increase in interlayer spacing, but also for the HER activity.</p><!><p>We have extensively (17 metals, 5 non-metals) studied substitutional doping of Mo and S in view of three complementary criteria: activity, feasibility and stability. The activity is assessed via hydrogen adsorption energies, with the promising systems being re-visited relying on grand-canonical DFT in order to account for the electrochemical potential. Furthermore, we have also quantified the improvement of the electrical conductivity of MoS 2 due to substitutional doping. The feasibility is addressed via two aspects: (a) relative stability of doping the edges vs. the basal plane and (b) the tendency of dopants to segregate (i.e. cluster). Finally, the stability of the proposed active sites is estimated by analyzing (a) the deformation induced by adsorbates, (b) release of H 2 S under HER conditions and (c) "blocking" of active sites by OH* adsorption. OH* adsorption can, however, also be seen as a modification of the active site that might still lead to HER, a possibility we have also assessed. Our computations identify 2Ti/S, 2Zr/S and 2Hf/S as the most promising subsitutional dopants in order to activate HER on the 2H basal plane of MoS 2 at overpotentials of about 0.2 V. For Mo/P, i.e. doping by phosphorus, the active site consists of OH*, but features a similarly small overpotential. Given our comprehensive analysis of substitutional doping and the few truly promising systems for activating the 2H basal plane of MoS 2 we expect to motivate experimental investigations to test these predictions.</p><!><p>This work was financially supported by Région Auvergne Rhône-Alpes through the project Pack Ambition Recherche 2018 MoSHi. We are very grateful to IFP Energies nouvelles for supporting the MoSHy project (N o 1801167601) and to the IFP Energies nouvelles team, Mona Marie Obadia and Quentin Cacciuttolo, for discussing the MoS 2 solids that should meet the predictions of molecular modeling in the case of hydrogen production by water electrolysis. The authors thank the SYSPROD project and AXELERA Pôle de Compétitivité</p>