Theoretical investigation of the non-metal sites of two-dimensional conjugated metal–organic frameworks based on benzenehexathiol for hydrogen evolution activity enhancement

Abstract

For electrocatalysts, the electrocatalytic activity of the non-metal sites is not negligible. We found that sulfur atoms should be the predominant active site for conjugated metal–organic frameworks (c-MOFs) based on benzenehexathiol (BHT) toward the hydrogen evolution reaction (HER). There is a “volcano”-shaped relationship between their HER activity and 3p band center of the sulfur active site. Interlayer interactions are also crucial in determining the HER activity of c-MOFs. Based on these findings, we proposed that Mo–BHT possesses excellent potential as an active HER catalyst.

---

10th Anniversary Statement

Journal of Materials Chemistry C has provided an excellent platform for high-quality research in electrically conducting metal–organic frameworks (MOFs). How to increase MOFs’ electrical conductivity has drawn much interest, especially regarding MOFs’ applications in energy-related fields. Many significant progresses about MOF-based electrochemical catalysts have been published here. Journal of Materials Chemistry C is a great journal that aims to publish innovative interdisciplinary studies, and that is why we are glad to submit this communication to this journal. We sincerely celebrate the 10th anniversary of Journal of Materials Chemistry C and wish it a bright future.

Electrocatalytic water splitting is an ecofriendly and sustainable way to produce hydrogen, which requires high-efficiency catalysts. Noble metals (Pt and its derivatives) are the most active electrocatalysts for the hydrogen evolution reaction (HER) in water splitting, but their high cost and scarcity greatly hinder their large-scale utilization.^1,2 Electrically conducting metal–organic frameworks (MOFs) are an emerging class of non-noble metal based HER catalysts.³ They can combine excellent electrical conductivity with the designable, tunable porous structure and chemical functionality of MOFs, providing great opportunities for the design of high-performance electrochemical catalysts.^4–11

Among electrically conducting MOFs, 2D conjugated metal–organic frameworks (c-MOFs) have superior charge transport properties, which are built by planar conjugated ligands like benzenehexathiol (BHT), hexahydroxytriphenylene (HHTP), etc. The electrical conductivity of Cu₃BHT can be up to 2500 S cm⁻¹.¹² In addition, 2D c-MOFs based on BHT (denoted as M–BHTs) not only have high electrical conductivity (Table S1, ESI†), but also contain the metal-dithiolene structural motif which is highly active for the HER in aqueous solution.^8,9,13 This makes it intriguing to explore the potential of M–BHTs as HER catalysts. Most M–BHTs have 2D lattices as shown in Fig. 1a. Early exploration of a M–BHT-based HER catalyst has been made by the Marinescu group.⁶ The HER performances of M–BHTs with similar structures but different metal centers have been explored, including Fe–BHT, Ni–BHT, and Co–BHT. Co–BHT exhibits a very low overpotential (185 mV) at a current density of 10 mA cm⁻², which is much better than those of Fe–BHT (473 mV) and Ni–BHT (331 mV).⁶ This shows that the HER performance of 2D c-MOFs is highly dependent on the metal center. However, systematic theoretical studies have found that the active site is not always located at the metal or sulfur sites of single-layer M–BHTs.¹⁴ In HER electrocatalysts like sulfide with metal–sulfur motifs, the sulfur site exhibits the highest active reactivity.^15–17 On the other hand, although the single-layer geometry of M–BHTs is commonly used in theoretical studies,^18–20 it probably leads to an incomplete understanding of catalyst activity since the information on the interlayer stacking structure is missing. In our previous study on the HER performance of Cu₃BHT, we found that Cu₃BHT nanoparticles, which expose more (100) surface, exhibit much better catalytic performance than the Cu₃BHT film and Cu₃BHT nanocrystals that expose the (001) surface.⁴ In order to explain this phenomenon, we used a three-layer model of Cu₃BHT to simulate the real catalyst, and we found that the most active sites on the extended (001) surface and edge site-rich (100) surface are different. These results inspire further insight into the source of the activity of c-MOFs in the HER, which will help us design a state-of-the-art non-noble metal HER catalyst.

Fig. 1 (a) Structure of a M–BHT (M represents metal atoms). The bilayer structures of (b) Co–BHT and (c) Mo–BHT (upper: top view; lower: side view). Blue, sulfur; grey, carbon; orange, cobalt; and red, molybdenum. Three kinds of S sites on the Mo–BHT surface are denoted as the S1, S2, and S3 sites.

In this study, the intrinsic activities and the electronic structures of the metal and non-metal sites of M–BHTs were comprehensively investigated. We also considered different metal centers (mainly for Co, Mo, and Ni and Fe) and interlayer interactions between M–BHT layers. The results of adsorption energetics in the HER process and electronic structure analysis showed that the sulfur site exhibits the highest inherent HER activity, making it the dominant active site. Instead of serving as an active site, the metal sites are more likely to affect the electronic structure of sulfur active sites, which explains why M–BHTs with different metal centers have significant differences in HER activity. We also found that the interatomic interactions between the metal site and the sulfur site in different M–BHT layers have marked influences on their HER activity. A relationship between the catalytic reactivity and the electronic structures of the sulfur active site was established. It provides valuable guidance for designing highly active M–BHTs catalysts from the perspective of tuning the surroundings of the non-metal site and considering interlayer interactions based on the multilayer M–BHT geometries. Finally, on the basis of these theoretical studies, we propose that Mo–BHT with a staggered bilayer structure should be a highly active HER catalyst.

In the HER process, the number of active sites and the intrinsic activity of each active site determine the overall performance of the electrocatalyst.²¹ The porosity and high surface area of M–BHTs catalysts (Fig. 1a) make them already expose a large number of active sites. Therefore, the inherent HER activity of M–BHTs should be dominated by the activity of each reaction site. Altering the metal center is an effective strategy to optimize the HER activity of M–BHTs. Co–BHT performs best among the M–BHTs, making it a vital benchmark (Table S1, ESI†). Mo is a potential metal center candidate to achieve such a goal because the Mo atom has more unoccupied d orbitals in the valence electronic configuration than the Co atom and thus has stronger electronic interactions with sulfur atoms. In addition, the doping of Mo in several non-noble metal based electrocatalysts significantly improves their HER activities,²² not to mention that MoS₂ is also a highly active HER catalyst.²³ Thus, the HER activity of Mo–BHT is worth exploring.

In the structure of a M–BHT, each 2D layer stacks with each other through van der Waals interactions. Therefore, the bulk structures of Co–BHT and Mo–BHT were considered to determine the more stable structures. The van der Waals interactions were described using PBE-D2 empirical corrections to account for long-range dispersion effects. After structural optimization (details can be found in Fig. S1 and Table S2 in the ESI† for the geometric parameters), we found that the Co–BHT layers are in an overlapping mode, while the Mo–BHT layers are in a staggered mode. Fig. 1b and c present the bilayer structures of Co–BHT and Mo–BHT, respectively. The stability and rationality of fully relaxed Co–BHT and Mo–BHT cells were demonstrated by the negative formation energies (Table S2 in the ESI†). To take into account the influence of adjacent layers on the HER activity, we adopted bilayer geometries for both Co–BHT and Mo–BHT (Fig. S2 in the ESI†). Specifically, the first layer was relaxed as the exposed (001) surface, and the second layer represented the Mo–BHT layer in the bulk phase. In Co–BHT, the Co, S, and C sites are identical for all the same kinds of sites. In Mo–BHT, the staggered stacking leads to strong interactions between S and Mo atoms belonging to different Mo–BHT layers, which has generated three kinds of S sites denoted as the S1, S2, and S3 sites (Fig. 1c).

The electronic structure of an electrocatalyst is important for understanding its HER performance.²⁴ For the electrocatalytic process of the HER involving proton–electron transfer, good conductivity of the catalyst is the foundation of its high catalytic performance.^4,25 As shown in the band structures of Co–BHT (Fig. 2a) and Mo–BHT (Fig. 2b), several energy bands cross the Fermi level, revealing that they are intrinsically metallic. The total density of states (TDOS; Fig. 2a and b) plots show that Mo–BHT exhibits a higher TDOS than Co–BHT at the Fermi level, revealing the better electrical conductivity.²⁶ The electron localization function (ELF) is a normalized parameter that ranges from 0 to 1, with higher values indicating a higher degree of electron localization. Here, we used the bilayer geometries of Co–BHT and Mo–BHT, which fully consider the Mo–S, S–C, and other possible interlayer interactions, to obtain accurate electron localization on the catalyst surfaces. The ELF plots of Co–BHT and Mo–BHT are presented in Fig. 2c and d, respectively. The ELF plots show that the strong Mo–S bond has caused electrons to delocalize over the whole conjugated backbone of Mo–BHT, which shows a higher degree of electron delocalization than Co–BHT. The differential charge density (Δρ) values were plotted between the integral M–BHT, the BHT unit, and the metal center through the equation: Δρ = ρ(MBHT) − ρ(C₆S₆) − ρ(M). In conclusion, there is effective electron transfer from the metal center to the S atom in the M–S bond in Mo–BHT (Fig. 2e and f). Therefore, the good conductivity of Mo–BHT and the effective interaction between the Mo centers and the organic linker guarantee its application in the HER as an active catalyst.²⁷

Fig. 2 Band structures and total density of states (TDOS) of (a) Co–BHT and (b) Mo–BHT. The dotted line represents the Fermi level. The electron localization function (ELF) of (c) Co–BHT and (d) Mo–BHT (the colorful bar is averagely segmented from 0 to 1). Differential charge density of (e) Co–BHT and (f) Mo–BHT with an isosurface value of 0.005 eÅ⁻³ (the yellow and cyan regions are representative of electron-rich and electron-defective regions, respectively).

The Gibbs free energy change was calculated to identify important issues in the HER process: the active site and reaction mechanism.²¹ There are two widely accepted reaction mechanisms for the HER: the Volmer–Heyrovsky pathway and the Volmer–Tafel pathway.²⁸ In these HER mechanisms, the HER process always starts with the initial state H⁺ + e⁻ to the final state 1/2H₂. The Gibbs free energy change of the intermediate state, an adsorbed H* (ΔG_H*, see the calculation details in the ESI†), has been generally accepted as a descriptor for the HER activity.^21,29 When ΔG_H* approaches zero, the HER activity is higher. Here, in order to determine the accurate active site of the M–BHT catalysts under consideration, all possible active sites on the surfaces of Co–BHT and Mo–BHT, including the metal site, the S site, and the C site, were investigated to find the “best” sites (Fig. S3 in the ESI†). The three-state Gibbs free energy diagrams for H* adsorption on Co–BHT and Mo–BHT at the equilibrium potential were calculated to estimate the intrinsic HER activity. As shown in Fig. 3a, the ΔG_H* values of the Co site and the C site on Co–BHT are larger than that on the S site (0.31 eV), indicating that the S site is the active site of Co–BHT in HER electrocatalysis. The dominant contribution of the S site on Co–BHT is similar to the MoS₂³⁰ and Cu₃BHT catalysts.⁴Fig. 3b summarizes the ΔG_H* values of all possible adsorption sites on Mo–BHT, namely the Mo site, the C site, and the S1, S2, and S3 sites, where these S sites are distinguished by different surroundings. Both Mo and C sites show high positive ΔG_H* values, revealing their poor contribution to HER performance. It can be seen that, depending on the bilayer structure of Mo–BHT, the ΔG_H* values of three S sites are different, which indicates the various activities of the HER. In the single-layer geometry of Mo–BHT, there is only one kind of S site, which should be used to determine the intrinsic activity. Thus, the single-layer geometry of M–BHTs results in an inability to determine the accurate active site, which highlights the significance of considering the interlayer interactions in theoretical calculations.

Fig. 3 Calculated Gibbs free energy diagrams for the HER on (a) Co–BHT and (b) Mo–BHT at the equilibrium potential. (c) Schematic diagram of the Volmer–Heyrovsky pathway and the Volmer–Tafel pathway of the HER mechanism on Mo–BHT, where a green ball represents the H atom (an isolated * denotes the adsorption site on the surface and H* represents an adsorbed H atom). Gibbs free energy diagrams for the HER following the Volmer–Heyrovsky pathway and the Volmer–Tafel pathway on (d) Co–BHT and (e) Mo–BHT.

The ΔG_H* of the S2 site, 0.29 eV, is the smallest value among all adsorption sites on Mo–BHT, which is smaller than that of Co–BHT. These results indicate that the S site is still the preferred reactive site for the HER on Co–BHT and Mo–BHT, regardless of the existence of various S sites induced by interlayer interactions. For the overall HER performance, although Co and Mo sites do not directly contribute significantly to the reactivity, they have made indirect promotion of the HER process on the S active sites through modulating the M–S interactions. Accordingly, the calculation results of ΔG_H* indicate that the S site contributes significantly to the overall activity of the HER, and Mo–BHT holds the potential to be a more active HER catalyst than Co–BHT. It is worth noting that the ΔG_H* of the S active site on Mo–BHT is smaller than that on Co–BHT. Therefore, Mo–BHT, supported by high conductivity, is believed to be capable of good HER catalytic activity. We expect the experimental evaluation of the HER performance of Mo–BHT to further confirm our hypothesis.

The reaction mechanism of the HER on Co–BHT and Mo–BHT catalysts was determined by calculating the free energy change along the Volmer–Tafel and Volmer–Heyrovsky reaction pathways. Take Mo–BHT as an example, the HER processes are presented in the schematic diagram in Fig. 3c. The calculation results show that the energy barrier of the Heyrovsky step is lower than that of the Tafel step in both Co–BHT (Fig. 3d) and Mo–BHT (Fig. 3e). For the whole HER process, the Volmer–Heyrovsky pathway is more thermodynamically favorable than the Volmer–Tafel pathway at the equilibrium potential. The Volmer step serves as the rate-determining step, which is in good agreement with the experimental results.^6,13 It can also be seen that the Gibbs free energy of the Heyrovsky reaction on Mo–BHT is lower than that on Co–BHT, which further highlights the potential of the Mo–BHT catalyst. Furthermore, the highest activity of the S active site and its direct interactions with reactive intermediates indicate that the effect of changing metal centers on the overall catalytic performance of M–BHTs is not from their own activity.

The highest HER activity of the S site and its leading role in the HER mechanism reveal the critical role of the S site in M–BHT catalysts. Thus, the electronic structures of the S site and metal site were further studied to elucidate their role in the origin of the electrocatalytic activity of M–BHT catalysts. Based on the modulation of the global electronic structure caused by changing the metal centers in M–BHTs, the projected density of states (PDOS) was calculated to distinguish the contributions of metal and S atoms. As shown in Fig. 4a and b, S 3p overlaps strongly with the d states of the metal atoms in both the bonding and antibonding regions. This property is conducive to enhancing charge delocalization in the MOF skeleton. As highlighted in the PDOS images, the electronic states of the Co 3d and S 3p orbitals of Co–BHT are larger than that of the C 2p orbital at the Fermi level, while in Mo–BHT, the Mo 4d orbitals make the most contribution in all orbitals. Correspondingly, by altering the metal centers in M–BHTs, the density of the electronic states of S sites varies in occupied states and at the Fermi level. More importantly, more electronic states of metal d orbitals at the Fermi level enable the adjacent metal atoms to regulate the electronic structures of the S active sites effectively. This was obviously reflected in the PDOS of three S sites on Mo–BHT, which will be discussed in detail later.

Fig. 4 The projected density of states (PDOS) of (a) Co–BHT and (b) Mo–BHT (the inset enlarges the PDOS near the Fermi level). (c) PDOS of the 3p orbital of the S1, S2, and S3 sites on Mo–BHT. (d) The relationship between −|ΔG_H*| and the 3p band center of S atoms in Mo–BHT, Co–BHT, Ni–BHT, and Fe–BHT.

The different environmental S sites induced by interlayer interactions provide a good platform for studying the relationship between HER activity and the electronic structure of S active sites. For the HER process on M–BHT, the higher the strength of the S–H* bond formed in the rate-determining Volmer reaction, the smaller the energy barrier in the critical step, and therefore the higher the intrinsic HER activity.³¹ The bond strength of S–H* depends on the match between the energy levels of the H 1s orbital and the S 3p orbital, where the reduction of the S 3p energy level can improve the matching degree. Mo–BHT just provides three kinds of S sites. The PDOS plots of the S1, S2, and S3 sites on Mo–BHT are shown in Fig. 4c. The distance between the S3 site and the Mo atom below in the second layer is 2.46 Å, indicating the formation of strong interactions between valence orbitals. The band center energies of the 3p electronic states of three S sites on Mo–BHT were calculated to provide a detailed comparison, which refers to the calculation method of the metal d-band center.³² The 3p band center energies of the S1, S2, and S3 sites on Mo–BHT are −2.70 eV, −2.73 eV, and −3.37 eV, respectively (Fig. 4c). As the S site attracts more electrons from adjacent metal atoms, the 3p band of the S site can be moved downward. The Bader charge analysis shows that the S1, S2, and S3 sites have −2.16 e⁻, −2.16 e⁻, and −2.20 e⁻, respectively. Traditionally speaking, the S3 site should possess the best HER activity. Yet, the HER activity does not always increase with the decrease of the 3p band of S sites over an expected wider energy range. To understand this anomaly, we also collected the ΔG_H* and 3p band center of the S sites of Co–BHT and the experimentally prepared Ni–BHT (AA and AB stacking modes) and Fe–BHT (see the details in Fig. S4 and Table S3 in the ESI†). It is found that the −|ΔG_H*| of M–BHTs and the corresponding 3p band center of S sites (S spin-down states) of Mo–BHT and Ni–BHT (Co–BHT and Fe–BHT) follow a “volcano”-shaped relationship (Fig. 4d). The equations of the fitted curves are as follows: y_L = 1.37x + 3.46 (the left curve) and y_R = −5.70 × −14.14 (the right curve), and y_L and y_R intersect at x = −2.49 eV. It can be seen that the S1, S2, and S3 sites are located on the left curve of the volcano, which clarifies the unusual HER activity of the S3 site of Mo–BHT. The results show that the electronic states of the S sites of M–BHTs can be affected by changing the metal atoms and establishing more interaction paths between the metal and S atoms. More importantly, the volcano-shaped relationship shows that the HER activity can reach the optimal ΔG_H* when the 3p band center of the S sites is modulated to −2.53 eV and −2.48 eV. This result points to a definite direction for climbing to the summit of the activity volcano and also provides a handy tool. The HER performance of novel M–BHTs can be evaluated and screened more efficiently by easy collection of the 3p band center of S sites.

In summary, the effect of the interlayer interactions between M–BHT layers on the HER activity origin of M–BHT catalysts was explored. Our calculation results reveal that the S atoms in the highly active Co–BHT and Mo–BHT are the catalytic active sites for the HER. By taking the stacking mode of M–BHT layers into consideration, differences in electronic states and the HER activity were found at S sites. In addition, a “volcano”-shaped relationship between the HER activity and the 3p band center of S active sites was established, which suggests that the optimal −|ΔG_H*| can be obtained when the 3p band center equals −2.53 eV and −2.48 eV. The various HER activities of S sites in Mo–BHT and Co–BHT also indicate that it is desirable to select favorable metal atoms and create more M–S interaction pathways to regulate the electronic states of S atoms to prepare excellent M–BHT catalysts. The pivotal roles of S atoms and metal atoms in the HER activity of M–BHTs were distinguished by theoretical studies. Our results also showed that considering the interlayer stacking mode of metal-dithiolene-based catalysts is required for an accurate evaluation of their electrocatalytic activities. We expect the favorable HER activity of Mo–BHT to be verified experimentally in the near future.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by the Scientific Research Foundation for High-level Talents of Anhui University of Science and Technology (2022yjrc102), the National Natural Science Foundation of China (Grant No. 52273170), the Strategic Priority Research Program of Chinese Academy of Sciences (XDPB13), and the Open Research Fund Program of Engineering Technology Research Center of Coal Resources Comprehensive Utilization, Anhui Province, Anhui University of Science and Technology (MTYJZX202201).

Notes and references

1. S. Lu and Z. Zhuang, Sci. China Mater., 2016, 59, 217–238.
2. A. Wang, L. Shen, M. Zhao, J. Wang, W. Zhou, W. Li, Y. Feng and H. Liu, J. Mater. Chem. C, 2019, 7, 8868–8873.
3. A. Nath, K. S. Asha and S. Mandal, Chem. – Eur. J., 2021, 27, 11482–11538.
4. X. Huang, H. Yao, Y. Cui, W. Hao, J. Zhu, W. Xu and D. Zhu, ACS Appl. Mater. Interfaces, 2017, 9, 40752–40759.
5. C. A. Downes and S. C. Marinescu, ACS Catal., 2017, 7, 8605–8612.
6. C. A. Downes, A. J. Clough, K. Chen, J. W. Yoo and S. C. Marinescu, ACS Appl. Mater. Interfaces, 2018, 10, 1719–1727.
7. R. Dong, Z. Zheng, D. C. Tranca, J. Zhang, N. Chandrasekhar, S. Liu, X. Zhuang, G. Seifert and X. Feng, Chem. – Eur. J., 2017, 23, 2255–2260.
8. R. Dong, M. Pfeffermann, H. Liang, Z. Zheng, X. Zhu, J. Zhang and X. Feng, Angew. Chem., Int. Ed., 2015, 54, 12058–12063.
9. C. A. Downes and S. C. Marinescu, J. Am. Chem. Soc., 2015, 137, 13740–13743.
10. A. J. Clough, J. W. Yoo, M. H. Mecklenburg and S. C. Marinescu, J. Am. Chem. Soc., 2015, 137, 118–121.
11. C. A. Downes and S. C. Marinescu, ACS Catal., 2017, 7, 848–854.
12. X. Huang, S. Zhang, L. Liu, L. Yu, G. Chen, W. Xu and D. Zhu, Angew. Chem., Int. Ed., 2018, 57, 146–150.
13. K. Chen, C. A. Downes, E. Schneider, J. D. Goodpaster and S. C. Marinescu, ACS Appl. Mater. Interfaces, 2021, 13, 16384–16395.
14. Y. Ji, H. Dong, C. Liu and Y. Li, Nanoscale, 2019, 11, 454–458.
15. T. Sun, J. Wang, X. Chi, Y. Lin, Z. Chen, X. Ling, C. Qiu, Y. Xu, L. Song, W. Chen and C. Su, ACS Catal., 2018, 8, 7585–7592.
16. J. Ye, Y. Zang, Q. Wang, Y. Zhang, D. Sun, L. Zhang, G. Wang, X. Zheng and J. Zhu, J. Energy Chem., 2021, 56, 283–289.
17. G. Li, T. Sun, H.-J. Niu, Y. Yan, T. Liu, S. Jiang, Q. Yang, W. Zhou and L. Guo, Adv. Funct. Mater., 2023, 2212514.
18. F. Li, X. Zhang, X. Liu and M. Zhao, ACS Appl. Mater. Interfaces, 2018, 10, 15012–15020.
19. J. Guan, T. Pal, K. Kamiya, N. Fukui, H. Maeda, T. Sato, H. Suzuki, O. Tomita, H. Nishihara, R. Abe and R. Sakamoto, ACS Catal., 2022, 12, 3881–3889.
20. Z. Ren, S. Wang, H. Zhang, B. Huang, Y. Dai and W. Wei, Phys. Chem. Chem. Phys., 2021, 23, 25239–25245.
21. Z. W. Seh, J. Kibsgaard, C. F. Dickens, I. Chorkendorff, J. K. Nørskov and T. F. Jaramillo, Science, 2017, 355, eaad4998.
22. J. Wang, T. Liao, Z. Wei, J. Sun, J. Guo and Z. Sun, Small Methods, 2021, 5, 2000988.
23. G. Ye, Y. Gong, J. Lin, B. Li, Y. He, S. T. Pantelides, W. Zhou, R. Vajtai and P. M. Ajayan, Nano Lett., 2016, 16, 1097–1103.
24. Y. Zheng, Y. Jiao, M. Jaroniec and S. Z. Qiao, Angew. Chem., Int. Ed., 2015, 54, 52–65.
25. L. Zhang, Z. Su, F. Jiang, L. Yang, J. Qian, Y. Zhou, W. Li and M. Hong, Nanoscale, 2014, 6, 6590–6602.
26. L. Sun, M. G. Campbell and M. Dincă, Angew. Chem., Int. Ed., 2016, 55, 3566–3579.
27. X. Wang, H. Niu, X. Wan, Z. Zhang, F. R. Wang and Y. Guo, J. Colloid Interface Sci., 2022, 624, 160–167.
28. B. E. Conway and B. V. Tilak, Electrochim. Acta, 2002, 47, 3571–3594.
29. J. K. Nørskov, T. Bligaard, A. Logadottir, J. R. Kitchin, J. G. Chen, S. Pandelov and U. Stimming, J. Electrochem. Soc., 2005, 152, J23–J26.
30. T. F. Jaramillo, K. P. Jørgensen, J. Bonde, J. H. Nielsen, S. Horch and I. Chorkendorff, Science, 2007, 317, 100–102.
31. Y. Jiao, Y. Zheng, K. Davey and S.-Z. Qiao, Nat. Energy, 2016, 1, 16130.
32. B. Hammer and J. K. Norskov, Nature, 1995, 376, 238–240.

---

Footnotes

† Electronic supplementary information (ESI) available: Computational details, the geometric parameters and stability of Co-BHT and Mo-BHT, calculation of the Gibbs free energy change, and additional figures and tables. See DOI: https://doi.org/10.1039/d3tc00238a

‡ These authors contributed equally.

---