New Insight on Hydrogen Evolution Reaction Activity of MoP2 from Theoretical Perspective.

We systematically investigated the hydrogen evolution reaction (HER) of six facets of MoP2 based on the periodic density functional theory (DFT). The calculated values of Gibbs free energy of hydrogen adsorption (ΔGH) indicated that the (111) facet has a good HER activity for a large range of hydrogen coverages. The zigzagged patterns before 75% hydrogen coverage suggest a facilitation among Mo1, P1 and Mo2 sites, which are attributed to repeat occupancy sites of H atoms. From ab initial atomistic thermodynamics analysis of hydrogen coverage, we gained that the most stable coverage of hydrogen is 18.75% at 1 atm H2 and 298 K. Finally, the doping effects on HER activity were investigated and found that catalytic performance can be improved by substituting P with an S or N atom, as well as substituting the Mo atom with an Fe atom, respectively. We hope this work can provide new insights on further understanding of HER for MoP2 and give instructions for the experimental design and synthesis of transition metal phosphides (TMPs)-based high-performance catalysts.

1. Introduction

class="Chemical">Hydrogen has been regarded as a class="Chemical">promising energy carrier and future sustainable high energy fuel, due to its class="Chemical">proclass="Chemical">perties of non-class="Chemical">pollution, energy class="Chemical">pan class="Chemical">conversion and gravimetric energy density [1,2]. The energetic efficiency, environmental sustainability and carbon-free alterative of hydrogen generations come from the interconversion between water and hydrogen. Electrolytic water-splitting [3] is one of the most reliable and effective methods to produce hydrogen in various hydrogen productions, which usually involves the hydrogen evolution reaction (HER). Therefore, it is very important to explore a stable, low cost, economically and environmentally friendly and high efficiency HER catalyst for realizing large-scale industrial electrocatalytic hydrogen production [4].

class="Chemical">Platinum and its alloys are class="Chemical">pan class="Chemical">considered to be the best active electrocatalyst to facilitate HER in acid condition because of its small overpotential and accelerating reaction rates [5], which can lead to more efficient water-splitting [6]. However, the high cost and low abundance of platinum limits its widespread application in practice [7,8]. Thus, it is highly necessary to explore earth-abundant and non-noble-metal HER electrocatalysts to substitute for platinum. Recently, many promising HER electrocatalysts have been unveiled, such as anatase [9], metal dioxides [10], carbides [11], selenides [12], sulfides [13], phosphides [14] and nitrides [15]. Among potential materials considered thus far, the family of transition metal phosphides (TMPs) has drawn a wide range of interests as a promising alternative catalyst for HER due to its high catalytic activities [16,17]. For instance, TMPs, including W-P [18,19], Fe-P [20,21], Ni-P [22,23] and Mo-P [24,25], are well known for their high catalytic activity and durability. Experimentally, many approaches were used to synthesize TMPs, such as the solvothermal method [26], hypophosphite thermal decomposition [27,28], high vacuum solid-state reaction method [28], metal-organic framework phosphatized method [29] and so on. Recently, researchers found TMPs with higher phosphorous component holds a higher catalytic activity for HER [29,30,31].

As a typical class="Chemical">TMP with higher class="Chemical">pan class="Chemical">phosphorous component, great stability and good electronical conductivity, has attracted tremendous attention from academia [32,33,34]. Experimentally, a facile two-step strategy was used to fabricate nanoparticle films on a metal Mo plate showing superior HER catalytic activity at all pH values [35]. Through a phosphidation reaction based on an nanosheet array on carbon cloth, an nanosheet was fabricated by Zhu et al. [36]. They found that the nanosheet delivers superior catalytic activity and stability with an overpotential as low as 58 mV to drive a catalytic current density of 10 mA cm−2. In addition, Wu et al. [28] synthesized nanoparticles with a high vacuum solid-state reaction method and carried out electrochemical tests. It was revealed nanoparticles show superior HER performance with small overpotentials of 38 mV and low Tafel slopes of 52 mV dev−1. More recently, Gao et al. [37] synthesized nanosheets using nanosheets as a precursor by high temperature solid state reaction with phosphorus under vacuum conditions. It was found that nanosheets exhibit very promising HER activity, considerably better than bulk and non-exfoliated . To our best knowledge, though there are many exciting experimental findings, insights into for the systematical mechanism of hydrogen absorptions are lacking, which is of critical importance for further enhancing the activities and stabilities of HER electrocatalysts.

Gibbs free energy of class="Chemical">hydrogen adsorclass="Chemical">pan class="Chemical">ption is usually used as a descriptor to predict HER activity of the electrocatalyst [38]. Nørskov et al. [39,40] reported that the maximum of HER activity is obtainable with the value of closing to zero. More negative values of mean the binding between hydrogen atom and facet is stronger, which is impedimentary for desorption of . On the contrary, more positive values of stand for the weaker blinding between hydrogen atom and facet. The step of proton/electron-transfer can be hindered with the positive value. Therefore, the best-performing and ideal electrocatalyst for HER should have a value of close to zero.

To figure out HER mechanisms of , the adsorclass="Chemical">ption structures and energetics of atomic class="Chemical">pan class="Chemical">hydrogen have been systematically calculated using periodic density functional theory (DFT). HER activities of the surface were predicted by the value of . Ab initio atomistic thermodynamics was employed to determine the most stable phase at 1 atm pressure and 298 K. Bader analysis [41] and different charge density [42] analysis were employed to further understand the interaction between Mo and P atoms. Finally, the doping effects on HER activity were investigated by substituting P with an S or N atom, as well as Mo atom with an Fe, Cu, Cr, Co, Mn or Ni atom, respectively. We hope our results can shed new light on further understanding of HER for and designing transition phosphide HER electrocatalysts.

2. Computational Methods

The spin-polarized periodic DFT calclass="Chemical">culations were class="Chemical">performed by the class="Chemical">plane-wave basis set with the class="Chemical">projector augmented class="Chemical">plan wave (PAW) method for treating class="Chemical">pan class="Chemical">core and valence electronics [43,44]. Perdew-Burke-Ernzerh (PBE) functional with generalized gradient approximation (GGA) [45] was carried to deal with the electronic exchange-correlation interaction. All calculations were implemented in the Vienna Ab initio Simulation Package (VASP) [46]. The plane-wave cutoff energy was tested using a well-known variable-controlling approach and set to 400 eV. The Gaussian smearing method was employed to represent the total energy with a smearing width of 0.02 eV. The low Miller-index facets were modeled by cleaved with the optimized structure of bulk MoP2, with a four-atom-layer 1 × 2 × 1 supercell for (100) facet and four-atom-layer 2 × 1 × 1 supercells for (111), (110), (101), (001) and (011) facets. To avoid the interaction among the adjacent images, 15 Å of vacuum along the perpendicular catalyst facets was added. The Monkhorst-Pack method with a centered k-point mesh of 3 × 3 × 1 was used to sample the Brillouin zone. The top half of the layers together with the added hydrogen atoms was relaxed during the geometry optimization in the whole calculations. During the structure optimization, the residual forces have converged to less than 0.025eV Å−1 and the total energy converged to less than 1.0×10−6 eV. To obtain a deep understanding of HER, the partial atomic charges were investigated using Bader charge analysis developed by Henkelman et al. [41]

The surface energies are determined as: where and are the total energy of the surface slab and the bulk, respectively. N is the number of formula unit in the slab and A is the facet area of the oclass="Chemical">ptimized slab. The difclass="Chemical">pan class="Chemical">ferential hydrogen adsorption energy is defined as: where and are the total energy of the system with n and (n-1) hydrogen atoms adsorbed on the surface. is the total energy of molecule in the gas phase. A negative value of suggests beneficial absorption. The differential Gibbs free energy of the adsorbed hydrogen atom is calculated by: here, is the adsorption energy of the hydrogen atom determined by Equation (2). and stand for the zero-point energy correction and entropy change between adsorbed hydrogen and hydrogen under standard conditions at the temperature of . can be obtained by ; here, is the entropy of an isolated hydrogen molecular in gas phase at standard condition. The value of is approximately equal to −0.2 eV. is defined as: where , and are the zero-point energy of n and (n-1) hydrogen atoms adsorbed on the surface, as well as a gas phase hydrogen molecular, respectively. The represents the zero-point energy of hydrogen molecule in gas phase. The zero-point energy can be calculated by: where h represents the Plank constant and is the vibrational frequency. The calculated value of the vibrational frequency is 4301 cm−1 for a single hydrogen molecular, which is consistent with the experimental finding of 4395 cm−1 [47]. The calculated values of vibrational frequencies are 1744.53 cm−1, 608.71 cm−1 and 548.19 cm−1 for hydrogen adsorbed on the (111) facet. Using Equation (4) and Equation (5), the values of 0.04 eV for was obtained for the (111) facet. Thus, the expression of can be represented as for the (111) facet.

To further identity the most stable phase, the Ab initio atomistic thermodynamics [48] was employed to determine at a specific temperature and pressure denoted by , which can be obtained by: where is the total energy of the system with the adsorbed class="Chemical">hydrogen atoms on the surface, is the total energy of the clean facet, is the total energy of moleclass="Chemical">pan class="Chemical">cule in gas phase and is the hydrogen chemical potential at a given temperature and pressure. Clearly, is a function of for a surface at a concrete hydrogen coverage. based on the specific condition can be determined by: here, is the hydrogen chemical potential at the standard pressure and is Boltzmann constant.

To investigate the kinetic mechanism of HER, the activation energies were calpan class="Chemical">culated based on the transition theory. The activation energy barrier () and reaction energy ( of the elemental reaction are obtained by the following formulas: where , and reclass="Chemical">present the total energy of initial state, transition state and final state, resclass="Chemical">pectively.

3. Results and Discussion

3.1. Electronic Structures of Bulk

Figure 1 shows the bulk , which has an orthorhombic structure with a space group of Cmc21 (36). A Mo atom is seven-class="Chemical">coordinated by six P atoms at ends of the class="Chemical">prism and one P atom at outside of the rectangular facet of the class="Chemical">prism, which shows class="Chemical">pan class="Species">five ring-pore structure and the zigzagging layered with a relatively large free volume. Table 1 lists the optimized lattice parameters that are in good agreement with the available experimental findings [28,49] and recent theoretical results [50]. The maximum error is less than 0.06 Å between the calculated and experimental findings.

Figure 1

Bulk . Mo: dark green; P: purple.

Table 1

The lattice parameters of bulk .

Bulk MoP_2.	a (Å)	b (Å)	c (Å)
Present	3.147	11.242	5.009
Calculation a	3.142	11.132	4.949
Experiment b	3.145 a	11.184	4.984

a. Reference [50] from density functional theory (DFT) with Perdew-Burke-Ernzerh (PBE) functional. b. Reference [49] from experimental data by X-Ray power method.

In surface science studies, priority will be given to low Miller-index facets due to their high durability. For orthorhombic , the low Miller-index facets can be cleaved to (100), (110), (111), (011), (101) and (001) examples. The oclass="Chemical">ptimized structures of these facets are class="Chemical">plotted in Figure S1 and class="Chemical">pan class="Chemical">corresponding surface energies are summarized in Table 2. Generally, the most stable facet has the lowest surface energy and can be synthesized in experiments [51,52,53,54]. From computing surface energy, it can be found the thermodynamic stability of the clean facets follows the order of (111) > (110) > (101) > (011) > (100) > (001). Thus, the most stable facet is (111), which agrees with the recent experimental finding where the (111) facet has been synthesized [54]. In addition, surface packing densities (SPD) are usually used to identify the stability. A larger value of SPD means that it is more stable. Table 2 also lists the calculated data of SPD for each facet. Likewise, the (111) facet is recognized as the most stable one. Therefore, the (111) facet is further examined.

Table 2

Calculated surface energies and surface packing density (SPD) of low Miller-index facets for .

MoP_2 Facet	(111)	(110)	(101)	(011)	(100)	(001)
Surface energy(meVÅ−2)	104.78	134.34	134.85	146.73	147.35	173.52
SPD (atom nm−2)	18	17	15	13	11	11

To further understand the electronic structure of the (111) facet, Figure 2 sketches the calclass="Chemical">culated total and class="Chemical">project electronic density of states (TDOS and PDOS, resclass="Chemical">pectively) data. Obviously, the facet shows metallic behavior with none-zero TDOS data at the class="Chemical">pan class="Chemical">Fermi level. It is advantageous for good conductivity which is beneficial for HER activity. The main contributions of TDOS are from Mo-d orbitals.

Figure 2

Project and total density of states of the (111) facet.

3.2. Hydrogen Adsorption

To investigate HER activity, all possible adsorclass="Chemical">ption sites on the (111) facet were carefully tested. It was revealed that there are six favorable adsorclass="Chemical">pan class="Chemical">ption sites located at Mo1, P1, bridge (Mo1-P3), Mo2, bridge (Mo1-P2) and P2, as shown in Figure S2. The absorption energies at each site were listed in Table S1, which follows the trend: Mo1 < P1 < bridge (Mo1-P3) < Mo2 < bridge (Mo1-P2) < P2. With the assumption that hydrogen atoms easily spread over the facet at the lowest energy sites for each type of hydrogen coverage, we depicted the global minimum of hydrogen adsorption for each type of hydrogen coverage on the (111) facet. Figure 3 shows the optimized structures with different hydrogen coverages (here, 100% coverage means adsorbed 16 H atoms on the (111) facet). It is usually thought hydrogen coverage should first occupy all sites with the lowest adsorption energies, i.e., all Mo1 sites on the (111) facet. Strangely, when the first H atom takes over one Mo1 site on the (111) facet, the second H atom do not go to the other Mo1 site. It occupies one P1 site and the third H atom occupies an adjacent Mo2 site, as shown in Figure 3a. Subsequently, occupancy sites of H atoms will repeat the procedure until reaching 75% coverage. For the last 75–100% coverage, additional H atoms prefer to adsorb on the bridge (Mo1-P3). Interestingly, there are no H atoms adsorbed on bridge (Mo1-P2) and P2 sites during the processes of hydrogen coverage. To further understand these phenomena, the adoption energies and charge density differences were investigated.

Figure 3

Optimized structures of the (111) facet with the H coverages: (a) 18.75%, (b) 37.5%, (c) 56.25%, (d) 75% and (e) 100%.

Figure 4 desclass="Chemical">cribes the calclass="Chemical">pan class="Chemical">culated values of and for H adsorption on the (111) facet with different hydrogen coverages. The absolute values of are less than 0.2 eV, which is a favorable value for HER. Therefore, it is indicated that the (111) facet can favor HER over a large hydrogen coverage range. The zigzagged patterns before 75% hydrogen coverage suggest a synergy among Mo1, P1 and Mo2 sites, which are attributed to repeating occupancy sites of H atoms.

Figure 4

Adsorption energy () and adsorption free energy ( ) as functions of hydrogen coverage on .

To obtain a deep understanding, herein, the charge density difclass="Chemical">ference induced by the adsorbed class="Chemical">pan class="Disease">H atoms on the (111) facet were analyzed. The charge density differences for H atoms adsorbed on the surface are defined by , where , and are the charge densities of the catalytic surface adsorbed H atoms system, the clean surface and H atoms, respectively.

Figure 5 shows the charge density difclass="Chemical">ferences for the system made uclass="Chemical">p of the (111) facet and the adsorbed class="Chemical">pan class="Disease">H atoms with H coverage of 1/16, 2/16 and 3/16. For 1/16 H coverage at Mo1 site, there are significant electron accumulation (yellow) and depletion (light blue) around the adsorbed H atom and the connected Mo atom, respectively. The election transfer mainly exists between H and Mo atoms, indicating a strong interaction between them. However, no electron depletion is found around the H atom, revealing that many electrons transferred from the Mo to the H atom. Interesting, a small amount of the electron depletion and accumulation are found around Mo1-P1 and P1-Mo2 bonds, respectively, which predict that the bond length of Mo1-P1 will increase, while that of P1-Mo2 will decrease compared with those in a clean facet, agreeing with the analysis of bond lengths from the geometry optimization in Table 3. Meanwhile, the small electron depletion of P atom at P1 site and accumulation of Mo atom at Mo2 site are in line with that from the Bader analysis in Table 3. Small amounts of charge density differences on the P1 site due to the occupation of the Mo1 site will affect the adsorption energy on the P1 site compared with that of only H coverage at the P1 site in Table 4, in which the adsorption energy will decrease. In other words, H adsorption at Mo1 site facilitates hydrogen adsorption at the P1 site, agreeing with a zigzag pattern with H coverage discussed above. For 2/16 H coverage, where the adsorbed H atoms are at both Mo1 and P1 sites, an accumulation (yellow) forms between H and P atoms due to a bonding generated between them. In addition, a depletion of the electron density outside the adsorbed H atom connected P atom occurs because of the rearrangement of electron toward the bond of H and P. When H coverage increase to 3/16, the phenomenon of electron accumulation and depletion around the H atom and connected Mo atom at Mo2 site are analogous to that of 1/16 H coverage.

Figure 5

The charge density differences for the adsorbed H atoms with H coverage of (a) 1/16, (b) 2/16 and (c) 3/16. Top and bottom panels represent overlooking and sidelooking, respectively. Color codes for Mo, P and H are dark green, purple and brown, respectively. Charge accumulation and depletion are plotted by the yellow and light blue regions with the isosurface value of 1.5×10−3 e/bohr3.

Table 3

Bond lengths and partial atomic charges from Bader analysis with H coverage of 1/16, 2/16 and 3/16.

	Clean facet	1/16 H atom	2/16 H atoms	3/16 H atoms
r_Mo_1−P_1 (Å)	2.313	2.330	2.380	2.394
r_p_1−Mo_2 (Å)	2.382	2.366	2.435	2.411
P1_Bader(e)	−0.360	−0.354	0.146	0.200
Mo1_Bader(e)	0.659	0.707	0.722	0.736
Mo2_Bader(e)	0.696	0.670	0.681	0.698

Table 4

Adsorption energies at the Mo1, P1 and Mo2 sites. Here, 1 means there is an adsorbed H atom, while 0 stands for 0 adsorbed H atoms.

Adsorption Energy	Mo1 site	P1 site	Mo2 site
ΔE_HeV	−0.106	0	0
	0	−0.038	0
	0	0	0.083
	1	−0.194	0
	1	1	−0.286

Overall, H adsorclass="Chemical">ption will lead to a rearrangement of electrons and change the bond length of class="Chemical">pan class="Chemical">Mo-P on the (111) facet shown in Table 3. Meanwhile, H adsorption at one site will facilitate hydrogen adsorption at the other site shown in Table 4.

Furthermore, to estimate the most stable class="Chemical">coverage of H on the (111) facet at 1 atm and 298K, ab initio atomistic thermodynamics was emclass="Chemical">ployed to identify the by the class="Chemical">pan class="Chemical">hydrogen chemical potential and pressure illustrated in Figure 6. Each line represents a given H coverage. Obviously, the most stable phase at 1 atm and 298 K is 18.75% H coverage, in which the active sites for HER is the Mo2 site. As shown in the Figure S3, the kinetic energy barrier of the Tafel pathway was considered by calculating the kinetic energy barrier profiles of the formation of Hydrogen, from which the HER activity of electrocatalysts can be known. The kinetic energy barrier profiles for the Tafel path at 18.75% H coverage is 0.75 eV, which is lower than the previous result of 1T- (0.85 eV) [41]. Based on the exiting learning, it can be learned that the lower energy barrier may facilitate the formation of the H-H bond.

Figure 6

at 298K with pressure and ranging from: (a) −1.4 to 0 eV and (b) −0.3 to −0.1 eV. Note that the dot line is used to clearly reflect the values of at 1 atm and 298 K. The intersection points show the values of for different coverages at 1 atm and 298K.

To understand the kinetic mechanism of HER at 18.75% H class="Chemical">coverage, the activation energies of the Taclass="Chemical">pan class="Chemical">fel reaction were calculated, in which two adjacent pre-bonded H atoms at top site of Mo2 and P1 combine to form an molecule. Figure S3 shows the kinetic energy barrier profile of HER on the (111) surface at 18.75% H coverage. Clearly, generation needs to experience two transition states. The first step needs to overcome a relatively high activation energy barrier of 0.75 eV. However, it drops to 0.18 eV for the second step. These low activation energy barriers revel that generation through Tafel reaction is favored on the (111) surface at 18.75% H coverage.

3.3. Doping on the (111) Surface

Doping has been extensively used to tune the physical and chemical properties of materials. It has been proven that doping on transition-mental is an efclass="Chemical">fective way, which influence the chemical activity of catalysts [55,56]. The class="Chemical">pan class="Chemical">hydrogen generation of catalytic activity of 2H can be promoted by the introduction of metal dopants such as Fe, Ni and Co [57,58,59]. Experiments found that N, S doping of MoP can enhance the HER activity [60,61]. To systematically recognize HER activity of the doping effects on HER activity for the (111) facet of have been done from a theoretical perspective.

Figure 7a illustrates the substitutional doing of the Mo site with a class="Chemical">Co, class="Chemical">pan class="Chemical">Cu, Mn, Ni, Fe or Cr atom, and doping of the P site with an S or N atom was evaluated. The Hubbard correction was considered for the metal dopant of Fe. Geometry optimizations of the doping system show that doping of the Mo site by using a Co, Fe, Mn, Ni or Cr atom, as well as doping of the P site by using a N or S atom could not induce visible structural change to the (111) facet. On the other hand, doping with a Cu atom induces slightly larger lattice distortion, which is 0.8% tensile strain, compared with the clear (111) facet. The calculated of the doped system was summarized in Figure 7b. Clearly, doping of the P site with an N or S atom and doping of the Mo site with an Fe atom give a value of closer to zero, which promote HER catalytic activity for the (111) facet of MoP2.

Figure 7

(a) Top view of the (111) facet with two doping atom sites. (b) Gibbs free energies of adsorption for the (111) facet and its doped facets.

4. Conclusions

In this work, HER catalytic activities for the (111) facet of were systematically investigated from theoretical perspective using periodic DFT calclass="Chemical">culations. The electronic structure, Gibbs free energy of class="Chemical">pan class="Chemical">hydrogen adsorption and ab initio thermodynamics of hydrogen coverage were calculated to reveal the HER characteristics of .

From the electronic structures of six low Miller-index facets, we found the thermodynamic stability of the clean facets follows the order of (111) > (110) > (101) > (011) > (100) > (001). As the most stable facet, the (111) facet has the metallic behavior with none having zero DOS at the class="Chemical">Fermi level. The values of indicated that the (111) facet can favor HER over a large class="Chemical">pan class="Chemical">hydrogen coverage range. The zigzagged patterns before 75% hydrogen coverage suggest a synergy among Mo1, P1 and Mo2 sites, which are attributed to repeat occupancy sites of H atoms. From charge density differences and Bader analysis, we found that H adsorption will lead to a rearrangement of the electron and change the bond length of Mo-P on the (111) facet. Meanwhile, H adsorption at one site will facilitate hydrogen adsorption at the other site. From the ab initio atomistic thermodynamics analysis, we obtained that the most stable phase at 1 atm and 298 K is 18.75% H coverage, in which the active sites for HER is Mo2 site. Finally, we found doping of the P site with an N or S atom and doping the Mo site with an Fe atom can promote HER catalytic activity for the (111) facet of .