CoP for hydrogen evolution: implications from hydrogen adsorption

Cobalt phosphide (CoP) is one of the most promising, earth-abundant electrocatalysts discovered to date for hydrogen evolution reaction (HER), yet the mechanism is not well understood. Since hydrogen adsorption is a key factor of HER activity, here we examine the adsorption of atomic hydrogen on the low-Miller-index surfaces of CoP, including (111), (110), (100), and (011), by using periodic density functional theory. From the calculated Gibbs free energy of adsorption, we predict that (111), (110), and (011) surfaces will have good catalytic activities for HER. From ab initio atomistic thermodynamics, we find that the stabilities of the surfaces at 1 atm H2 and 300 K follow the trend of (111) 4 (100) B (110) c (011). On the most stable (111) surface, both Co bridge sites and P top sites are found to be able to adsorb hydrogen with a close-to-zero free energy change and the synergy of proximal Co and P atoms on the surface results in a better HER activity. Our work provides important insights into CoP’s excellent HER activity and a basis for further mechanistic understanding of HER on CoP and other transition-metal phosphides.


Introduction
Water splitting using renewable energy can provide a sustainable supply of fuel for future societies with hydrogen as a key energy carrier. 1 To date, platinum remains the most efficient electrocatalyst for hydrogen evolution reaction (HER). However, the low natural abundance and high cost of platinum hamper its wide use at the industrial scale. Thus, it is highly desirable to develop efficient, low-cost and earth-abundant electrocatalysts for HER and an enormous amount of research efforts have been devoted to it over the past decade.
Two types of possible pathways have been proposed for the mechanism of HER in acid media: the Volmer-Heyrovsky vs. the Volmer-Tafel mechanism. The Volmer reaction refers to the initial adsorption of protons from the acid solution to form adsorbed H (H + + e À -H ad ). 47 In the Volmer-Heyrovsky mechanism, a solvated proton from the water layer reacts with one adsorbed surface hydrogen to form H 2 (H ad + H + + e À -H 2 ), 48 while in the Volmer-Tafel mechanism, two adsorbed surface hydrogens next to each other react to form an H 2 molecule (H ad + H ad -H 2 ). 49 Thus, good HER electrocatalysts should be able to attract protons from the solution, while still can desorb H 2 . Nørskov and coworkers have correlated experimental exchange currents with calculated hydrogen adsorption energies from density functional theory (DFT), where maximum activity is obtainable when Gibbs free energy of hydrogen adsorption (DG H ) is zero. 50 This is mainly because lower DG H results in a very strong binding of the atomic hydrogen and hinders desorption of H 2 , while higher DG H prevents the binding of atomic hydrogen on the catalyst surfaces. Both of them will slow the reaction. Therefore, DG H is considered as a good descriptor of HER activity, and good HER electrocatalysts are expected to have close-to-zero DG H .
To understand HER on CoP in particular and to provide insights and guidelines for designing transition metal phosphide HER electrocatalysts in general, herein we have studied in detail the adsorption structures and energetics of atomic hydrogen on several low-Miller-index surfaces of CoP from first principles DFT. We use the calculated DG H to predict the HER activities of the surfaces. We further employ ab initio atomistic thermodynamics to determine the most stable and active surface of CoP at 1 atm H 2 pressure and 300 K.

Computational
Spin-polarized DFT calculations were performed by using the Vienna ab initio simulation package (VASP). 51 The ion-electron interaction was described with the projector augmented wave (PAW) method. 52 Electron exchange-correlation was represented by the functional of Perdew, Burke and Ernzerhof (PBE) of generalized gradient approximation (GGA). 53 A cutoff energy of 400 eV was used for the plane-wave basis set. (111), (110), (100), and (011) surfaces were examined. About 10 Å thick slabs in (2 Â 2) lateral cells with 15 Å of vacuum along the z-direction were used to model the adsorbate-surface systems for (111), (110), (100), and (011) surfaces; the Brillouin zone was sampled by (3 Â 3 Â 1) Monkhorst-Pack k-point mesh. Both the layer thickness and k-point mesh were tested to achieve a convergence of hydrogen adsorption energy within 0.04 eV. The top half of the slab was allowed to relax together with the adsorbed H atoms and the convergence threshold for structural optimization was set to be 0.025 eV Å À1 in force. Partial atomic charges were obtained using Bader charge analysis as implemented by Henkelman and co-workers. 54 We also used a different method 55 to derive the atomic charges based on the electrostatic potential, called the ESP charges.
Surface energies was determined by using where E slab is the total energy of the surface slab, E bulk is the total energy of the bulk unit cell (eV per CoP), N is the number of unit formula in the slab, and A is the surface area. The differential hydrogen adsorption energy DE H was calculated by where E(CoP + nH) and E[CoP + (n À 1)H] represent the total energy of the CoP system with n and n À 1 adsorbed hydrogen atoms on the surface, respectively, and E(H 2 ) represents the total energy of a gas phase H 2 molecule. A negative value of DE H suggests favorable absorption. The differential Gibbs free energy of adsorption DG H was obtained by where DE ZPE is the difference in zero point energy between the adsorbed H and H in the gas phase H 2 molecule, and DS H is the entropy difference between the adsorbed H and H 2 in the gas phase at standard conditions. To locate the most stable configuration for each coverage of H on a surface, we allowed the H atoms to first find the lowest energy sites and then fill the next lowest-energy sites and continue until the desired coverage was reached. In other words, we focused on finding the global minimum of hydrogen adsorption, with an assumption that the diffusion of H atoms on the surface would be facile and they always go to the lowest-energy sites. Ab initio atomistic thermodynamics 56 was used to identify the most stable phase of each surface at 1 atm H 2 and 300 K. The Gibbs free energy of adsorption at a specific temperature and pressure DG ad (T,p) was calculated by where E total (N H ) is the total energy of the CoP system with N H hydrogen atoms adsorbed on the surface, E total (0) is the total energy of the clean CoP surface, and E total H 2 is the total energy of a gas phase H 2 molecule. Dm H (T,p) is the chemical potential of hydrogen at different temperatures and pressures. Eqn (4) allows one to plot the Gibbs free energy of adsorption as a function of the hydrogen chemical potential for a specific surface at a specific H coverage. At any given hydrogen chemical potential, the stability of these H-covered surfaces can then be compared.

Dm H (T,p) is related to specific (T,p)-conditions by
where Dm H (T,p 0 ) is the hydrogen chemical potential at the standard pressure which can be looked up from the thermodynamic tables. 57

Results and discussion
We start with the bulk structure of CoP, followed by the clean surfaces including (111), (110), (100), and (011). We then examine and compare hydrogen adsorption on these surfaces.

Bulk CoP
Bulk CoP has a B31 or MnP-type structure, with a symmetry group of Pnma. In this type of structure, each metal atom is surrounded by six phosphorous atoms situated in a distorted octahedral configuration (Fig. 1a), while each phosphorous atom is surrounded by six metal atoms forming a highly distorted triangular prism (Fig. 1b). The experimental 58 and calculated lattice parameters of bulk CoP are listed in Table 1. One can see that they are in good agreement.

Clean surfaces
Low Miller-index surfaces are usually considered first in surface science studies due to their higher stabilities than higher index surfaces. We consider here four such surfaces: (111), (110), (100), and (011). Fig. 2 shows the structures of these clean CoP surfaces and Table 2 shows their calculated surface energies. One can see that the stabilities of the clean surfaces in vacuum follow the trend of (011) 4 (110) 4 (111) 4 (100). In general, surfaces with higher atom packing densities are expected to have higher stabilities. Here, we find that the surface packing densities of (011), (110), (111), and (100) are 25, 24, 24, and 22 atom per nm 2 , respectively, which indeed confirms the highest stability of the (011) surface and the lowest stability of the (100) surface. To learn the electronic structures of these surfaces, their electronic density of states (DOS) have been calculated. As shown in Fig. 3, they are all metallic and have non-zero DOS at the Fermi level. Being the most stable one, the (011) surface also has the lowest DOS at the Fermi level.
Hydrogen adsorption on the CoP (111) (111) is that DG H shows a zigzag pattern with H coverage from 25% to 75%, which closely correlates with the adsorption sites, as discussed next. Fig. 5 shows the optimized structures of CoP(111) with different hydrogen coverages. We find four types of stable hydrogen adsorption sites on CoP(111). For the first 25% hydrogen atoms, they are strongly adsorbed at cobalt bridge sites. During 25-75% hydrogen coverage, there is an alternative adsorption between cobalt bridge sites and phosphorous top sites, leading to a zigzagged DG H with a small fluctuation. More important, DG H is close-to-zero during this range, so both cobalt bridge sites and phosphorous top sites could be the active sites for HER on the (111) surface. For the 75% to 100% coverage, the extra hydrogen atoms are adsorbed at cobalt top sites, leading to more than one hydrogen atoms on some cobalt atoms.
The ''zigzag'' pattern between 25% and 75% hydrogen coverage suggests a synergy between proximal cobalt and phosphorous sites in adsorbing H atoms, conductive to HER. To understand this behavior, Table 3 shows differential hydrogen adsorption energies (DE H ) at the phosphorous top site at 0% and 31% hydrogen coverages. One can see that when there is no hydrogen adsorbed on Co (0%), the interaction of H with P is weak but when there are hydrogens already adsorbed on Co (31%), the H-P interaction becomes stronger. In other words, hydrogen adsorption at Co facilitates hydrogen adsorption at P. We found that this is due to the weakening of the Co-P bond induced by hydrogen adsorption on Co, as evidenced by the increasing Co-P bond length, which results in less negative charge on P (both Bader 54 and ESP 55 charges show the same trend) and stronger H-P covalent bonding. This synergy between Co and P for adsorbing H, as shown in Scheme 1, could be a key to CoP's high HER activity.
To evaluate the most stable coverage of H on CoP(111) at the experimental HER conditions (1 atm H 2 and 300 K), we used ab initio atomistic thermodynamics to determine DG ad (T,p) as a function of hydrogen chemical potential (Dm H ) or pressure at T = 300 K. Fig. 6 shows the calculated DG ad (T,p) vs. Dm H relationship: each line represents a given H coverage. One can Fig. 1 Bulk CoP: (a) each Co atom is surrounded by six P atoms situated in a distorted octahedral configuration; (b) each P atom is surrounded by six Co atoms at the corners of a highly distorted triangular prism. Co, green; P, red.  View Article Online see that at extremely low hydrogen pressure, clean surface is the most stable, while at very high hydrogen pressure, 100% hydrogen coverage is the most stable (Fig. 6a). At ambient pressures, Fig. 6b shows that 75% hydrogen coverage is the most stable with DG ad = À19.53 meV Å À2 at 1 atm. At this coverage, both Co and P sites adsorb H.
Hydrogen adsorption on the CoP(110) surface For a favorable adsorption (DG H o 0), we found that DG H is close to zero when hydrogen coverage is in the range of 25-50% (DG H B À0.10 eV). When hydrogen coverage is above 50%, DG H is slightly positive but also close to zero. Fig. 7 suggests that the CoP(110) surface can also be good for HER activity for a large range of H coverage. Fig. S1 (ESI †) shows the optimized structures of CoP(110) at different hydrogen coverages. From 0% to 100% coverage, hydrogen atoms occupy cobalt bridge   sites, phosphorous top sites, cobalt top sites, and another kind of cobalt top sites successively. Fig. 8  This might be the reason why we do not find DG H B 0 on CoP(100). Fig. 9b plots DG ad as a function of Dm H for CoP(100) and we find that 33% hydrogen coverage is the most stable at 1 atm and 300 K with DG ad = À13.56 meV Å À2 .
Hydrogen adsorption on the CoP(011) surface Fig. 10a shows DG H and DE H for CoP(011) at different hydrogen coverages. One can see that the adsorption of H on CoP(011) is in general much weaker than on the other three surfaces. On the other hand, this leads to a quite large range of H coverages from 0 to 50% where DG H B 0, indicating that CoP(011) may have good activity for HER. Fig. S3 in ESI † shows the optimized structures of CoP(011) surface at different hydrogen coverages. We find two types of stable hydrogen adsorption sites on CoP(011). The first 50% coverage of hydrogen atoms prefer phosphorous top sites with DG H B 0, so phosphorous top sites could be the active sites for HER on the (011) surface. This is in agreement with a recent study. 33 The last 50% coverage of Table 3 Differential hydrogen adsorption energies (D E H ) at the P site, Co-P bond length (r Co-P ), and partial atomic charges on P at 0% and 31% hydrogen coverages from Bader analysis (P Bader ) and electrostatic potential (P ESP )   hydrogen atoms occupy cobalt bridge sites. Fig. 10b    sites and P top sites are found to be able to adsorb hydrogen with a close-to-zero free energy change, indicating that the synergy of proximal Co and P atoms on the surface can result in a better HER activity. Third, hydrogen coverage, pressure, and temperature all have effects on the stabilities of the surfaces. We find that the stabilities of the clean surfaces in vacuum are different from the stabilities at 1 atm H 2 and 300 K. The stabilities of the clean surfaces in vacuum follow the trend of (011) 4 (110) 4 (111) 4 (100), while the stabilities of the surfaces at 1 atm H 2 and 300 K follow the trend of (111) 4 (100) B (110) c (011). This suggests that if one purposefully prepares the catalyst with the (111) facet, then the catalyst is likely to have a longer lifetime given its high stability during the reaction conditions, while the (011) surface might transform to other more stable surfaces under reaction conditions.
Although the present work focuses on the thermodynamics of HER on different CoP surfaces from computed H adsorption free energy, our results also have interesting implications for the reaction rates. For example, 75% coverage of hydrogen on CoP(111) would be very helpful for the Volmer-Tafel mechanism, as higher surface coverage would increase the chances for two adsorbed hydrogens next to each other to react and form an H 2 molecule. We plan to study the detailed mechanism of HER on CoP surfaces to examine the effect of coverage on reaction kinetics, as recently done for HER on MoS 2 . 59

Conclusions
To shed light on the activity of hydrogen evolution reaction (HER) on CoP, we have studied the adsorption of atomic hydrogen on the low-Miller-index surfaces of CoP, including (111), (110), (100), and (011), by using periodic DFT. From the calculated DG H , we predict that (111), (110), and (011) surfaces will have good catalytic activities for HER. From ab initio atomistic thermodynamics, we find that the stabilities of the surfaces at 1 atm H 2 and 300 K follow the trend of (111) 4 (100) B (110) c (011). So combining DG H and stability, we conclude that the (111) surface of CoP is the most promising facet for high and long-lived HER activity. The key feature of hydrogen adsorption on the (111) is that both Co bridge sites and P top sites are able to adsorb hydrogen with a close-to-zero free energy change and that there is a synergy of proximal Co and P atoms on the CoP(111) surface in adsorbing hydrogen. These insights open a door for further mechanistic understanding of HER on CoP and other phosphides.