H₂ holds great promise for meeting the ever-growing energy needs of humankind, and decreasing greenhouse gas emissions; however, significant efforts are required to develop H₂ production technologies that would allow to achieve its full potential [1, 2]. Water electrolysis is considered to be the ideal method for the production of H₂, because it can be powered using intermittent renewable energy sources, particularly wind and solar energy, and it forms a zero-carbon-emission energy cycle when associated with the fuel cell technology [3, 4]. Although promising, owing to its high overpotential, the critical half reaction, known as the hydrogen evolution reaction (HER), requires the use of noble metal catalysts [5-7]. However, the use of noble metal catalysts, such as Pt and Pt alloys, greatly increases the cost of the process and thus limits its widespread application [8, 9].

To develop low-cost and highly active HER catalysts, scientists have focused their efforts on transition metals which could be used as HER catalysts [9еC16]. Transition metal phosphides, such as CoP, WP, and CuP, have attracted the attention of scholars owing to their remarkable catalytic activity [13еC16]. Among these materials, CoP presented high catalytic activity in both acidic and alkaline media, and Tian et al. [13] designed a type of self-supported CoP nanowire arrays, which performed well in the pH range of 0 to 14. Du et al. [14] designed CoP nanotubes via template-assisted synthesis with the Faradaic yield of approximately 100%. Moreover, it was revealed that when CoP was decorated with other transition metals, it outperformed pure CoP [15, 16]. Density functional theory (DFT) calculations indicated that Zn doping optimized the thermal-neutral hydrogen adsorption of CoP, and doped CoP exhibited better HER activity than pristine CoP [15]. Zhang performed DFT calculations and revealed that the thermal-neutral free energy of Al-doped CoP was higher than that pure CoP, and Al-doped CoP reached the HER current density of 10 mAcm⁻² at the low potential of 23 mV [16]. Although several studies have been conducted on the preparation of doped CoP catalysts, their catalytic mechanism remains unclear. That hinders to obtain the material design principle, and hence, catalyst preparation becomes a "trial and errorА" process.

The structure of WP is similar to that of CoP; moreover, WP presents HER catalytic activity. However, theoretical calculations indicated that the HER catalytic performance of WP would be worse than that of CoP, because the adsorption of H by WP would be too strong. According to our calculation on the (101) facet of CoP, H adsorption on a fraction of activity sites, for which H desorption was limited, was too strong. Therefore, we speculated that if the surface of CoP were doped with W, the strong interactions between W and H would weaken the adsorption of H on the Co sites owing to the poisoning effect, which could promote the HER on CoP [17]. In addition, the structure similarities between CoP and WP indicated that W doping was more feasible than doping using other transition metals, and therefore, W-doped CoP could be a promising non-precious metal HER catalyst. In our preliminary study, we synthesized W-doped CoP and confirmed that its HER catalytic activity was superior to that of pristine CoP [18]. In this study, DFT was further utilized to investigate the HER mechanism of W-doped CoP catalysts based on the following questions: (1) How does W influence the HER catalytic activity of CoP? (2) What is the suitable W content to improve the HER catalytic activity of CoP? (3) Which are the suitable sites that W should occupy? This study provides a different method for obtaining doped CoP and could be helpful for designing more effective doping approaches.

The generalized gradient approximation with Perdew, Burke, and Ernzerhof of the DMol3 model of the Material Studio 2017 package was used to describe the electron exchange and correlation [19, 20]. DFT semi-core Pseudopots was used for core treatment [21]. The basic set was DNP 4.4 (double-numeric quality basic set with polarization functions) and the TS method was used for DFT-D correction [22]. The spin was unrestricted during calculations; moreover, 5 × 8 × 4 and 3 × 3 × 1 Monkhorst-Pack grids were used to sample the bulk CoP (WP) and slab models. The slab models consisted of six layers which were built using 1 × 2 supercells. The geometry convergence tolerances for energy change, max force, and max displacement were 1 × 10⁻⁵ Ha, 0.002 Ha/Å, and 0.005 Å, respectively. The smearing of 0.008 Ha was used. Thermodynamic calculations revealed that the structures were stable.

The adsorption energy (E_ad) , which was used to estimate the kinetics of the HER, was calculated from the free energy using Equation (1) [23, 24]:

(1)

Moreover, the free energy (G) was calculated as follows [24]:

(2)

where E is the total energy of system, T is the temperature, S is the entropy, and ZPE is the zero-point energy. ZPE was calculated as ∑(hv_i/2), where h is the Planck constant and v_i is the vibrational frequency [25]. The total free energy of the H⁺ ions and e⁻ was equal to that of 1/2H₂ at pH 0 when the potential was 0.

The HER is a two steps process that consist of the adsorption and desorption of H. Taking the d band center into account during the analysis of the HER is very useful [23, 26]. The center of the d band could be calculated using the density of states as follows:

where x is the energy and p(x) is the density of state. According to the d band theory, if the d band center of the adsorption site is low, the binding strength between the active sites and H atoms is also low, because more electrons are loaded into the antibonding orbitals, and therefore the desorption of H atoms is favored [23, 26].

The (101) facets of CoP and WP were used to build the slab model in this study; this was determined to be the main exposed surface, and was considered to be the active facet for the catalysis of the HER [15, 16, 27]. According to our preliminary study, only two types of Co sites were present on the surface of CoP. Compared to surface or subsurface doping, when W was doped in the body of CoP, its structure energy increased, which indicated that near-surface W doping was preferred. This tendency could be further explained by the atomic radius of W being larger than that of Co, and therefore, it would be more difficult to insert W into the body of CoP from a thermodynamics perspective. Furthermore, the effect of W doping would be weakened as the W sites would be far from the surface. When W is doped at different sites of the same CoP slab, the structure energies are very similar (< 0.5 eV), and therefore, all CoP structures that were surface or subsurface doped with different amounts of W were considered in our study.

The basic model of CoP was built according to the data reported by Selte, Birkeland, and Kjekshus (ICSD #624584) [28]. After geometry optimization, the lattice parameters were determined to be a = 5.080 Å, b = 3.280 Å, and c = 5.590 Å (Fig. 1(a)). Two types of Co atoms (marked as Co1 and Co2) and two types of P atoms (marked as P1 and P2) were present on the surface of CoP (Fig. 2(a)). The deformation charge density of CoP is illustrated in Fig. 1(b), and the results were consistent with the Mulliken charge analysis (Co1 lost 0.436 e⁻ and Co2 lost 0.463 e⁻; see Table 1). The structures of WP and CoP are similar. They were modeled using the data reported by Guerin, Sergent, and Prigent (ICSD #42056) [29], and the results are presented in Fig. 1(c). Two types of W atoms (marked as W1 and W2) and two types of P atoms (marked as P1 and P2) were present on the surface of WP (Fig. 2(b)). The surface of bulk CoP was marked as type I. The surfaces of CoP that contained 8.4 wt%, 16.8 wt%, and 33.6 wt% W were marked as types II, III, and IV, respectively. Surface or subsurface W doping did not change the structure of CoP significantly.

In the slab model denoted as top-view in Fig. 2(a), two types of Co atoms were present in each layer; the first type, Co1, was in line with two P atoms, and the second type, Co2, was located at the edge. Moreover, the two different types of P atoms (P1 and P2) in each layer were classified according to their side-view (Fig. 2(a)). For the WP model, the W and P atoms were marked using the same method. The d band centers of Co1, Co2, W1, and W2 were calculated to be −2.212, −2.235, −2.741, and −3.002 eV, respectively (Fig. 2(c) and (d)).

There are twelve types of possible adsorption site on the type I surface, namely four top sites (T_Co1, T_Co2, T_P1, and T_P2) and eight bridge sites (B_Co1Co1, B_Co2Co2, B_P1P1, B_P2P2, B_Co1P1, B_Co1P2, B_Co2P1, and B_Co2P2). According to the structure of CoP, the distances between the P atoms and other atoms were too long to form bridge bonds during the adsorption of H atoms [30]. After relaxation, five feasible sites remained, namely three top sites (T_Co1, T_Co2, and T_P2) and two bridge sites (B_Co1 and B_Co2) as presented in Fig. 3 and Table 3. The free energy of the H adsorption process (△G) is the key parameter of the HER, and reflects its feasibility [18, 19]. First, the closer △G is to zero, the easier it is for the H atoms to be desorbed, which is beneficial for the HER. If △G is too high or too low, it is difficult for the H atoms to be adsorbed or desorbed, which is detrimental for the HER. As summarized in Table 3, the △G values indicated that the most stable site was B_Co2 (−0.36153 eV). In addition, T_Co2 (0.09336 eV) and T_P2 (0.09798 eV) were determined to be more beneficial for the reduction of H because the corresponding △G was close to zero, which would promote the desorption of H, and thus would facilitate the transfer of protons [23]. In contrast, B_Co1 and T_Co1 presented more negative △G values during the adsorption process, which indicated that it was more difficult for H to be desorbed from the B_Co1 and T_Co1 sites. The low △G values of the WP sites in Table 3 indicated that it was difficult for the H atoms to be desorbed from most sites, and that would limit the HER.

According to the △G values of B_Co1 (−0.36153 eV) and T_Co1 (−0.21875 eV) the Co1 atoms attract H atoms, while the Co2 atoms repel H atoms (△G > 0 eV). The difference between T_Co1 and T_Co2 could be explained using the partial density of state (PDOS, Fig. 2(c)) because the d band center of the Co1 atom (−2.0094 eV) was higher than that of the Co2 atom (−2.0813 eV), and therefore the adsorption of the Co1 atoms toward H atoms was stronger than that of the Co2 atoms. The large d band difference between the Co1 and Co2 atoms might be caused by the different geometric structures around them. Therefore, the Co1 and Co2 atoms were marked as valid and invalid sites, respectively.

Figure 4 depicts the decorated W-doped CoP models, where one surface Co atom is replaced by a W atom (8.4 wt% mass, type II); these atoms were denoted W1-CoP and W2-CoP (when W replaced Co1 and Co2 atoms, respectively) according to their decoration sites. The d band center of the surface Co atoms decreased after W doping, as illustrated in Fig. 4(c) and (d). According to the charge distribution data presented in Table 4, the decrease in the d band center could be explained by the stronger electronegativity of W than that of Co, which could limit the gain of electrons on P (Table 4). As a consequence, Co would lose less electrons and its d band center would decrease, which would be beneficial to the desorption of H atoms.

The △G values of the feasible sites were further calculated and the results are listed in Table 5. The △G values of most W1-CoP sites were lower than those of CoP, which indicated that the H adsorption capability of those sites was enhanced. In detail, the △G of T_Co1 decreased to −0.25639 eV, the △G of T_Co2 was maintained at approximately +0.093 eV, and the △G values of T_P2 and T_P2ʹ became closer to zero. The bridge sites disappeared after W doping. For Co1, the disappearance of the bridge site was the result of the large difference in H adsorption between Co1 (−0.25639 eV) and W1 (−0.45860 eV). Moreover, the T_W1 sites could become inactive toward the HER because it was difficult for H atoms to be desorbed from these sites. The balance of the bridge site of Co2 was broken. The T_P2 and T_P2ʹ sites were more suitable for the HER because their △G values became closer to zero. B_P sites did not exist because the distances between P atoms were too large for them to form bonds. In summary, when one Co1 site was occupied by one W atom, the H desorption at the other sites would be intensified, and therefore, the HER would be promoted on the W1-CoP model.

When the Co2 atoms were replaced with W atoms, the model was denoted W2-CoP. The strongest site of CoP was B_Co1. According to the △G values of the W2-CoP sites listed in Table 5, △G of B_Co1Co1 was −0.23098 eV, and the main adsorption sites presented weaker adsorption than CoP did (−0.36153 eV). Moreover, the △G of the T_P2 site of W2-CoP was lower △G (−0.13805 eV) than that of the T_P2 site of pristine CoP. Similar to the situation of W1-CoP, T_W2 was determined to be the favorite site for H atoms and the △G of the T_W2 site of W2-CoP (−0.35801 eV) was lower than that of the B_Co1 site of CoP (−0.36153 eV). This indicated that the desorption of H on W2-CoP was intensified and the HER would be promoted via a different mechanism.

After partly replacing Co with W in CoP, the changes in △G could be explained by the changes in the d band center of the adsorption sites. First, as illustrated in Fig. 4(c) and (d), the insertion of W into CoP lowered the d band center of the Co atoms at the surface, and led to the increase in the △G values of the Co sites (the T_Co2 and B_Co2Co2 sites of W1-CoP and the B_Co1Co1 site of W2-CoP) [23, 31]. Second, H atoms prefer W over Co according to the △G values of the feasible sites of WP and CoP (Table 3), and thus, H atoms would be adsorbed on the T_W sites and the △G of the other sites adjacent to the W sites (the T_Co1, T_P2, and T_P2ʹ sites of W1-CoP and the T_P2 site of W2-CoP) would decrease when the Co1 or Co2 atoms were replaced by W atoms. Lastly, the above-mentioned effect mostly depended on the decoration site: when W was inserted into an invalid site (Co2), the adsorption of H atoms at the T_W site of W2-CoP was weakened (−0.35801 eV); in contrast, the attraction for H atoms of the W sites was activated (−0.45860 eV) when the valid site (Co1) was replaced by W.

To further confirm the effect of W on the catalytic process, two W atoms were used to occupy the Co1 and Co2 sites (type бз) in three models (2W1CoP, 2W2CoP, and W1W2-CoP), as illustrated in Fig. 5. The W atoms at the Co1 and Co2 sites were denoted as W1 and W2, respectively. Moreover, the PDOSs of the d band are depicted in Fig. 5(d) and (e). It was determined that the d band center of the surface Co atoms decreased when two W atoms were inserted.

According to Table 6, when two Co1 atoms were displaced, the △G values of T_Co2, T_W1 and B_W1W1 were +0.12345, −0.53316, and −1.04221 eV, respectively. When two Co2 atoms were displaced, the △G values of T_Co1, B_Co1Co1, T_W2, and B_W2W2 were −0.27104, −0.34603, −0.17492, and −0.17756 eV, respectively. When a Co1 site and a Co2 site were occupied by two W atoms, T_Co2 exhibited strong repulsion toward H atoms (+0.25059 eV). The feasible sites were B_Co1W1, T_W2, T_W2, T_P2 and T_P2ʹ (−0.47180, −0.48511, −0.47113, −0.07691, and +0.08067 eV, respectively). It would be more difficult for H atoms to be desorbed from the 2W1-CoP and W1W2-CoP sites, because the △G of the main metal site indicated that its adsorption capacity was too strong, which would be harmful for the HER. In contrast, the 2W2-CoP sites would be more suitable for the HER, as their △G would be closer to zero on similar site. The △G of P sites changed very little when the W doping content was 16.8 wt%.

As the prediction of valid site and invalid site, the W1 performs strong adsorption of H that even approaches the adsorption effect of W sites in WP (Table 3). However, when two W atoms were located at invalid sites, the adsorption capacity of W would decrease significantly. When one type of site (valid or invalid) was displaced, the desorption on the other type of site would be enhanced, and that was attributed to the decrease in the d band center. After doping, the P1 sites remained inert and adsorption on the P2 sites was slightly enhanced.

To evaluate the properties of CoP with high W doping content (33.6 wt%), all Co sites on surface were occupied (W4-CoP, type Ⅳ, Fig. 6). The △G values of the feasible surface sites are presented in Table 7. Among these sites, T_P2 presented the lowest △G (−1.25173 eV); moreover, the △G of B_W1W1 (−0.91225 eV) was lower than that of B_W2W2 (−0.30399 eV) and the △G of T_W1 (−0.42734 eV) was lower than that of T_W2 (−0.25291 eV). These results were consistent with the above-mentioned tendency that valid sites presented stronger adsorption than the invalid ones. Thus, the adsorption capacity of the CoP crystals could be significantly enhanced by decorating CoP with high W amounts.

Considering that exposure to W was detrimental to the desorption step of the HER, we used W atoms to replace Co atoms on the subsurface of CoP crystals, as illustrated by the six models in Fig. 7. The Co sites located in line with the P sites in the side view are denoted as Co1 and the other Co sites are denoted as Co2. The replacement of one Co1 atom and one Co2 atom with one W atom each is depicted in Fig. 7(a) and (b), respectively; moreover, the replacement of two Co1 and two Co2 atoms with two W atoms each is presented in Fig. 7(c) and (d), respectively. Lastly, the replacement of one Co1 site and one Co2 site by W atoms and the substitution of all Co sites on the surface of CoP with W atoms are illustrated in Fig. 7(e) and (f), respectively.

Table 8 lists the △G values of the feasible adsorption sites of the six above-mentioned models. When a W atom was inserted in a subsurface Co1 site (denoted as W1), the adsorption on most sites (T_Co2, T_Co2ʹ, and T_P2) was enhanced. The HER catalytic activity of the surface P1 site adjacent to W1 was enhanced, and a bridge bond was formed between the Co1 site and a H atom. The disappearance of the T_Co1 and B_Co2Co2ʹ sites was accompanied by the appearance of the B_Co1P1 site, and that implied the weakening of the adsorption capacity of the metal sites, which was in agreement with the decrease in their d band center. The slightly improved adsorption on the B_Co1Co1 site might be caused by the geometrical effect attributed to the insertion of W atoms. According to the △G values of the model 2 sites listed in Table 8, the situation was similar to that of model 1. With the exception of the △G of B_Co1Co1ʹ, those of most sites (T_Co1, T_Co2, T_Co1ʹ, B_Co2Co2, T_P2, and T_P2ʹ) increased. The data for model 3 suggested that when two subsurface Co1 sites were occupied, compared to model 1, the desorption on most feasible surface sites would be further enhanced. For model 4, two W2 atoms located at Co2 sites could improve the adsorption of H on the surface Co1 sites. When two W atoms replaced subsurface Co1 and Co2 atoms (model 5), most top metal sites (T_Co2, T_Co1ʹА», and T_Co2ʹ) disappeared and the bridge sites between Co and P (B_Co2P1, B_Co1ʹP1, and B_Co2ʹP2ʹ) appeared. The catalytic activities of the T_Co1, B_Co2P1, and B_Co1ʹP1 sites were superior to those of the other sites. For model 6, four Co sites were replaced by W atoms, and the adsorption capacities of the metal sites was intermediate between those of models 3 and 4.

In summary, subsurface W doping at the dopant content of 8.4 wt% could improve the HER catalytic activity of the surface by enhancing the desorption of H atoms. Further increasing the dopant content to 16.8 wt%, the two W atom located at the Co1 sites presented the best catalytic activity of all sites, owing to the △G values of all feasible sites being closer to zero. When the W content increased to 33.6 wt%, the HER catalytic activity of the W-doped CoP was the worst because the number of feasible sites decreased and the △G values of most sites were far from zero.

The best adsorption sites of the six above-mentioned models were always located near valid sites (Co1), and their △G values were −0.05710, −0.02539, −0.07488, −0.01247, and −0.20796 eV. This implied that the valid sites were the main HER sites of CoP.

Although numerous doping experiments were conducted to improve the activity of CoP catalysts, theoretical studies on the doping effect and catalysis mechanism are still lacking. In this work, a first-principles study was performed, and the results revealed that two types of metal active sites were present on the surface of the (101) facet of CoP, where the Co atoms at the valid sites (Co1) attracted H atoms and the Co atoms at the invalid sites (Co2) holds exclusion. After doping, it is determined that W played a double-sided role in the modification of CoP. On one hand, W doping could reduce the d band center of Co at the surface, which would be beneficial for the desorption of H atoms. Conversely, surface adsorption was intensified partly because H atoms preferred to bond with W owing to its strong adsorption capacity. Moreover, because the W atoms located at valid sites presented strong attraction toward H atoms and the adsorption was weakened at the invalid sites, it was suggested that W atoms performed similarly to Co at the valid and invalid sites. Furthermore, the low content of W dopant could improve the HER catalytic activity of CoP by improving the desorption capacity of the crystal surface and increasing the number of feasible sites. However, the HER catalytic activity would decrease as the content of doped W increased because of the strong adsorption between W and H atoms. Lastly, when using W doping to reduce the d band center of Co, inserting W atoms into the subsurface of CoP was an effective method for weakening the detrimental effects of W. Therefore, the HER catalytic activity of surface-doped CoP with 8.4 wt% W was much higher than that of pristine CoP, and the HER catalytic activity of CoP could be further improved via subsurface doping with 16.8 wt% W. The insertion of W atoms at sites located far from the valid sites was more effective, and therefore, selective doping could increase the HER catalytic activity of CoP catalysts.

The authors declare no competing financial interest.