Toward a deeper understanding of H–Ni₃C interactions: rule-based insights

Abstract

Ni₃C exhibits high catalytic activity for hydrogen evolution and hydrogenation reactions; however, the electronic structure of the Ni₃C surface has not been sufficiently investigated to clarify its interaction mechanism with H. In this study, we theoretically elucidate the fundamental rules governing the interaction between H atoms and Ni₃C surfaces: (1) the formal charge and valence of C are −2 and 2, respectively, (2) bonding Ni to C has a formal charge of +2, and (3) the formal valence of Ni²⁺ is 2. Resonance structures based on these rules enable an a priori discussion of the geometries, electronic structures, and interactions of Ni₃C surfaces with H atoms. Specifically, reconstruction of the Ni₃C(113) surface, correlation between Ni–H bonding strength and the number of Ni–C interactions, and heterolytic dissociation of hydrogen have been derived, which are important for understanding Ni₃C catalysis. These insights cannot be obtained from the d-band centre theory, nor can they be completely predicted using state-of-the-art machine learning interatomic potentials. This study clearly demonstrates the importance of understanding the degree of metal–carbon covalency for accurately interpreting the catalytic activity of nickel carbides.

---

1. Introduction

The effective utilisation of hydrogen is important for achieving sustainable development goals (SDGs). Hydrogen fuels produce no carbon dioxide emissions, and catalytic hydrogenation of organic compounds is a fundamental reaction in the synthesis of chemicals and pharmaceuticals. Therefore, catalysts with low energy and high selectivity should be developed for water electrolysis and hydrogenation.

Interactions between metal surfaces and adsorbates are typically interpreted using the d-band centre theory,^1–3 which predicts that hydrogen bonding weakens as the d-band centre shifts downward and vice versa. Since the d-band centre is determined by the lattice constant and valence electron count of the metal, catalytic activity for hydrogen evolution and hydrogenation can be optimised by alloying, forming a core–shell structure, or nanostructure control. However, this approach is similar to perturbation, and catalyst development based on the d-band centre theory approaches its limits. Consequently, new strategies for enhancing catalytic performance are required.

In recent years, nanomaterials based on metal carbides have attracted considerable attention for heterogeneous catalysis. In contrast to metal oxides, metal carbides do not have bandgaps and exhibit metal-like properties in general.^4,5 Additionally, metal carbides exhibit high charge polarisation and induce larger changes in the density of valence electrons on metals than in alloy systems. Hence, metal carbides can be used as electrode materials for water electrolysis and nanoporous metal carbides and nanoparticles exhibit excellent catalytic activity.^6–20 In particular, nanoparticles of nickel carbide (Ni₃C) show high catalytic activity for gas- and liquid-phase reactions, such as the dry reforming of methane and hydrogenation for organic synthesis in the liquid phase, as well as electrode reactions, such as hydrogen evolution.^9–18 Mizugaki et al. experimentally revealed that the activation barrier for hydrogenation of polar functional groups by Ni₃C/Al₂O₃ is a quarter that of Ni/Al₂O₃ (Ni₃C: 14.9 kJ mol⁻¹, Ni: 59.8 kJ mol⁻¹).¹⁷ Thus, carbidation of metals is a new promising way of improving catalytic activity, leading to more efficient use of hydrogen.

However, the origin of the catalytic activity of metal carbides remains under debate: can their high activity for the hydrogen be fully explained by d-band centre theory? does the p-block element, C, play a role in the catalytic mechanism? Density functional theory (DFT) calculations of carbides of group 4–9 elements, such as Ti, V, Co, Fe, Zr, Mo, and W, have been performed to answer these questions.^21–43 Nevertheless, the universal theory for understanding the catalytic activities of metal carbides has not been suggested. Furthermore, the isolation, synthesis, and especially the catalytic functions of group 10 metal (Ni, Pd, Pt) carbides have lagged behind, resulting in limited fundamental understanding of the electronic structures of Ni₃C despite its high catalytic activity.^9–18 Therefore, we assume that systematising the interaction between Ni₃C and hydrogen is important for systematising the catalytic activity of metal carbides. Temperature programmed desorption of H₂ (H₂-TPD) is effective for experimentally analysing hydrogen interactions. However, in H₂-TPD studies of Ni and Ni₃C nanoparticles, significant dependences on their supports are observed.^17,44,45 Even on Al₂O₃ supports, which are generally assumed to be difficult to hydrogen spillover, a support effect is confirmed for Ni₃C catalysts.¹⁷ Although some studies have investigated the catalytic activities of Ni₃C by DFT calculations,^15,16,18,46 to the best of our knowledge, no work has provided a detailed analyses of Ni₃C–H interactions, based on electronic structure analysis, as well as a clear theoretical framework for understanding these interactions.

Understanding the unique catalytic activity of Ni₃C requires comprehension of Ni₃C–H interactions based on its electronic structure. Adsorbed H on Ni₃C surfaces is an important intermediate in the catalytic reactions of Ni₃C nanoparticles, such as hydrogen evolution, methane reforming (where H₂ is a product), and liquid-phase hydrogenation. To suggest a universal theory for understanding Ni₃C–H interactions, the present study has investigated the electronic structure of Ni₃C in detail using DFT calculations. A rule to describe the electronic structure of Ni₃C has been derived from the calculated results, and based on the theory of resonance structures incorporating this rule, we have shown that the properties of Ni₃C including its interaction with H can be explained. We have compared Ni₃C with other representative metal carbides and discuss whether the rules established for Ni₃C could be the first step towards understanding the catalytic activity of metal carbides. Furthermore, the effectiveness of machine learning interatomic potential (MLIP)^47–49 for predicting the Ni₃C–H interaction based on the rule and the usefulness of MLIP in achieving theory-driven and data-driven studies on metal carbides are discussed. The discussion clarifies the issues and potential of the state-of-the-art MLIP.

2. Computational procedure

The PBE functional was used as the exchange–correlation functional in the DFT analysis.⁵⁰ A plane-wave basis set was used for the calculations, and the core region was treated using the projector-augmented wavefunction (PAW).^51,52 The numbers of valence electrons in Ni, C, and H were 10, 4, and 1. The energy cut-offs for the wave function and augmented charges were set to 600 and 2400 eV, respectively. Electronic states were optimised using self-consistent field calculations with an energy threshold of 10⁻⁷ eV. The geometry optimisation threshold was set to 10⁻² eV Å⁻¹. For the bulk model, both lattice parameters and atomic positions were optimised, whereas for the slab model, only the atomic positions were optimised. The geometry of the slab models was optimised with the bottom layer fixed to mimic the bulk structure (single-sided slab models were used), unless otherwise noted. DFT calculations were performed using the VASP.^53–56

For the bulk model of Ni₃C, the calculations were conducted using a k-point sampling of a 12 × 12 × 4 Monkhorst–Pack mesh;⁵⁷ the dependence on k-point sampling is shown in Fig. S1 (SI). Based on the dependence on bulk Ni₃C, the k-point sampling of the Ni₃C slab model was set to 4 × 6 × 1. For calculations involving molecules, a 3 × 3 × 3 nm³ supercell was used and only the Γ point was sampled. A 12 × 12 × 12 Monkhorst–Pack mesh was used in the calculations for Ni. Calculations were also carried out with dispersion force correction (DFT-D3(BJ))^58,59 and on-site Coulomb correction (DFT+U);⁶⁰ however, the effects were found to be negligible (Fig. S1 and S2). The present study discusses the adsorption structure of hydrogen atoms; hence, the influence of dispersion forces is negligible. Conversely, dispersion force correction will be crucial for investigating the adsorption structures of aromatic molecules. The negligible effect of dispersion correction on the structure and electronic states of Ni₃C is an important finding for future works that investigate adsorptions of large molecules. The dependence on on-site Coulomb parameters may be significant when comparing other metal carbides; this is presented in Section 3.4. Therefore, unless otherwise noted, this study uses the calculated results without applying on-site Coulomb correction and dispersion force correction for discussion.

The Ni₃C surface was modelled as a slab with the (113) plane. The (113) plane has been adopted in all previous studies performing DFT calculations on Ni₃C.^15,16,18,46 Furthermore, XRD measurements of the Ni₃C nanoparticle catalyst show that the (113) plane has the highest intensity.^15–18 Although there is no direct experimental evidence that (113) is the most exposed or active surface, considering the comparison with previous studies and the XRD peak, the (113) plane was adopted as the model surface in this study. The Ni₃C(113) slab models have a five Ni-layer thickness and a vacuum region of >17 Å. The bottom layers of slabs were fixed during optimisation for mimicking the bulk structure.

The number of valence electrons per atom was estimated using the Bader charge analysis algorithm^61–64 and the atomic charges were defined as the difference between the number of electrons in the valence region of the PAW method. Projected crystal orbital Hamilton population (pCOHP) analysis was performed using the LOBSTER package.^65–67 Using the projected atomic orbitals by LOBSTER, Mulliken and Löwdin population analyses were performed.⁶⁸ The dependence of populations on the methods is summarised in Table S1 (SI): there is no qualitative difference. The electron density distributions were visualised using VESTA.⁶⁹

In recent years, energy and structure predictions using MLIPs have been increasingly used in high-throughput data-driven research as MLIP calculations are 100–1000 times faster than DFT, while maintaining comparable accuracy.^47–49 The application of MLIP in materials chemistry research is in its infancy, and benchmarking studies are needed to assess their accuracy. MLIP calculations were therefore performed using the CHGNET potential,⁴⁸ a model suitable for electrode materials, and the preferred potential (PFP), a universal potential for surface reactions;⁴⁹ however, the correlation between the DFT and CHGNET calculations was poor (Fig. S3, SI). The correlation between the PFP and DFT data is presented and discussed in Section 3.5. PFP version 5.0 was used and calculations were performed in the CRYSTAL_U0 mode.

3. Results and discussion

3.1. Electronic structure of bulk Ni₃C

Fig. 1 shows the bulk Ni₃C unit obtained after geometry optimisation. The optimised lattice constants (a = b = 4.595 Å, c = 12.995 Å) are in good agreement with experimental values^70,71 and previous DFT calculations^42,71,72 (see Table S2). The C atoms occupy the octahedral sites of Ni with a coordination number of 6. In bulk Ni₃C, the Ni atoms are bridged by C atoms, and the resulting topology is shown in Fig. 1(c). The Ni–Ni distance in bulk Ni₃C is 2.66 Å, which is wider than that in bulk Ni metal (2.49 Å).

Fig. 1 (a) Perspective and (b) top views of the optimised structure of bulk Ni₃C. (c) Topological view of Ni₃C. Green and black balls represent Ni and C atoms, respectively.

The atomic charges of the Ni and C atoms in Ni₃C were compared with those in Ni, NiO, CH₄, and CO (Table 1). The Ni in Ni₃C is oxidised, although its formal valence is lower than that in NiO (+2). The absolute value of the atomic charge of C in Ni₃C is closer to that in CO than to the zero valence of C in CH₄, indicating a formal valence of −2 for C in Ni₃C. These results suggest that Ni₃C exhibits ionic character.

Table 1 Bader atomic charges of Ni and C atoms in Ni₃C, Ni metal, NiO, CH₄, and CO

Atom	Compound	Charge [e^−]	Atom	Compound	Charge [e^−]
Ni	Ni3C	+0.42	C	Ni3C	−1.28
Ni	Ni	0.00	C	CH4	−0.04
Ni	NiO	+1.04	C	CO	+1.83

---

Fig. 2 shows the calculated density of states (DOS) of bulk Ni₃C. Since Table 1 indicates its ionic nature, we examined the effect of the on-site Coulomb parameter on the DOS (Fig. S2); however, no significant changes were observed that would impact the discussions in this study. The total DOS (Fig. 2(a)) reveals a continuous band of Ni₃C across the Fermi energy (E_F) level, indicating metallic behaviour, consistent with experimental observations.^9,71 No spin polarisation was found in the DOS distribution of the α (major spin) and β (minor spin) electrons. Therefore, unlike Ni, bulk Ni₃C is not expected to exhibit spontaneous magnetisation, thus aligning with previous reports.^9,72Fig. 2(b) presents the projected DOS of Ni and C atoms. The band at the Fermi level originated from Ni states, and no C states were identified around the E_F. C states can be identified below −4.9 eV. A clear correlation between the C and Ni orbitals is observed in the energy range of −4.9 to −7.5 eV, whereas little Ni–C orbital correlation is observed from the peak at −12 eV. These Ni and C states were projected onto the orbital angular momentum component (Fig. 2(c) and (d)). The PDOS of Ni (Fig. 2(c)) confirms that the band at the Fermi level is the d component of Ni. In the Ni–C correlation band, Ni had a small s component in addition to d components. In bulk Ni, bands originating from the 4s orbitals were observed, whereas in Ni₃C, no s components were identified in the bands around the Fermi level, except for the Ni–C correlation band. From the correlation with C, it is confirmed that the s component of Ni contributes significantly in the energy range from −7.5 to −5.8 eV, while the d component of Ni dominates in the range from −4.9 to −5.8 eV. The split in the Ni–C correlation band was due to the difference in the orbital angular momentum component of Ni. The PDOS of C (Fig. 2(d)) reveals that the Ni–C correlation band (−4.9 to −7.5 eV) has a p orbital contribution from the C atom, and the independent band (−11.5 to −13.5 eV) originates from the s orbitals of the C atom; thus, the states of the s and p components of the C atoms are completely separated, and it is concluded that no s–p hybridised orbitals such as sp, sp², and sp³ are constructed in Ni₃C. Fig. S3 (SI) shows the DOS projected onto the magnetic angular momentum components. No significant splitting was identified in the Ni d orbitals (d_xy, d_yz, d_zx, d_x²−y², d_z²) or in the C p orbitals (p_x, p_y, p_z).

Fig. 2 Density of states (DOS) of bulk Ni₃C. (a) Total DOS, (b) projected DOS of Ni and C atoms, (c) projected DOS of the s, p, and d components of Ni, and (d) projected DOS of the s and p components of C. E_F is the Fermi energy.

We now discuss the electronic structure of Ni₃C as derived from the calculated results. First, the ionisation states of the Ni and C atoms are considered. The valence-electron configuration of the Ni atom is 3d⁸4s², and the 4s electrons are transferred to the 2p orbitals of C atoms upon ionisation. As illustrated in Fig. 3(a), this results in a formal valence of +2 for Ni and −2 for C. C²⁻ has two unpaired electrons that can be used for covalent bonding, whereas Ni²⁺ has two unpaired electrons in the 3d orbital. Hence, the formal valences of C²⁻ and Ni²⁺ are 2, and they form a covalent bond, as shown in the PDOS results. In addition, even when Ni does not covalently interact with C, the d orbitals of Ni can participate in the Ni–Ni metal bonding (Fig. 3(a)). Therefore, the bands originating from the Ni d orbital remained on Ni₃C, and the d band of Ni₃C was located independently in the Ni–C correlation band. This d band in Ni₃C across E_F (E − E_F > −4 eV in Fig. 2(b)) has similar properties to the d band of bulk Ni, whereas the DOS on E_F is reduced by carbidation because part of the d orbital participates in the bonding interactions with the C atoms. Based on the Stoner model,⁷³ the spontaneous magnetisation of bulk Ni is due to exchange interactions between different electronic spins, and spin polarisation is induced to resolve the repulsion at E_F as the DOS increases. The decrease in the DOS at E_F in Ni₃C resulted in the disappearance of spontaneous magnetisation, as indicated by the total DOS of Ni₃C (Fig. 2(a)). In addition, the Ni–Ni bonds in Ni₃C were weaker than those in the Ni metal because of the decrease in the DOS of the Ni d band, contributing to Ni–Ni metallic bonding, and the Ni–Ni bond lengths in Ni₃C are longer than those in the Ni metal (Fig. 1). Thus, the formal Ni²⁺/C²⁻ interaction summarised in Fig. 3(a) adequately illustrates the geometry and band structures of Ni₃C. Specifically, the three frontier bands in Ni₃C are the (1) Ni–C correlation band, primarily composed of Ni 4s and C 2p orbitals (−5.8 < E − E_F < −7.5 eV, corresponding to the ionic interaction between Ni and C), (2) Ni–C correlation band, primarily composed of Ni 3d and C 2p orbitals (−4.9 < E − E_F < −5.8 eV, corresponding to the covalent interaction between Ni and C), and (3) the d band of Ni (E − E_F > −4 eV, corresponding to the metallic interaction).

Fig. 3 (a) Schematic of the electronic configuration of C²⁻ and Ni²⁺ in bulk Ni₃C. (b) Resonance structure of Ni₃C.

The C atoms in Ni₃C are hexacoordinated, and the symmetrically equivalent C atoms are bridged by Ni atoms. The resonance structures of the Ni₃C bond (Fig. 3(b)) are derived from its electronic configuration (Fig. 3(a)) and topology (Fig. 1(c)). Ni, which bonds to C (green in Fig. 3(b)), bridges two C atoms, as in –C–Ni–C–, when the formal valence charges of Ni and C are set to +2 and −2, respectively, forming an infinite covalent bond, which corresponds to the Ni–C correlation band shown in Fig. 2(b). Because all the Ni atoms in Ni₃C are equivalent, the resonance can be represented as shown in Fig. 3(b), resulting in a formal Ni valence of +2/3. This explains the results of the Bader analysis, in which the valence of Ni in Ni₃C is lower than that in NiO. As shown in Fig. 3, the Ni in Ni₃C exhibits metallic band character. Additionally, due to the contribution from the resonance structure, the formal valence of Ni is expected to be +2/3 or less. In fact, the XPS spectrum of Ni in Ni₃C approximates that of Ni foil, and no distinct Ni²⁺ peak, such as that observed in NiO, is detected.¹⁶ The resonance structure (Fig. 3(b)) derived from the electronic configuration (Fig. 3(a)) is the key point of this study. The covalent interaction between Ni and C is important for understanding the properties of Ni₃C. In the following sections, we verify the theory that the properties of Ni₃C surfaces can be described by the resonance structures depicted with the following rule:

(1) The formal charge of C is −2 and its formal valence is 2.

(2) The Ni bonded to C²⁻ is assumed to be Ni²⁺.

(3) The valence of Ni²⁺ is 2.

3.2. Electronic and geometric structures of Ni₃C(113)

The (113) plane, with the highest intensity peak in the XRD spectrum of the Ni₃C nanoparticles, was used for discussion. Previous studies have assumed that the (113) plane will be the most exposed and catalytic-active surface of the Ni₃C nanoparticles, and investigated the interaction with H atoms.^15,18,46 There are two possible terminations in Ni₃C(113) (Fig. 4(a) and (b)), and researchers have argued that the C-terminated surface is more stable.⁴⁶ DFT calculations were performed for dual-sided slab models, and the surface energies of the C-terminated and Ni-terminated Ni₃C(113) were estimated (Fig. S5 in the SI). The calculated results are consistent with those reported in a previous study.⁴⁶ As shown by the DOS of the C-terminated and Ni-terminated surfaces (Fig. 5), where surface states are found between the Ni d band and the Ni–C correlation band. The PDOS results for the Ni₃C models are summarised in the SI (Fig. S6–S8).

Fig. 4 (a) C-terminated Ni₃C(113), (b) Ni-terminated Ni₃C(113), and (c) reconstructed Ni₃C(113). Green and black balls represent Ni and C atoms, respectively.

Fig. 5 Total DOS of C-terminated Ni₃C(113) (a) and Ni-terminated Ni₃C(113) (b). Arrows represent the surface states. E_F is the Fermi energy.

The PDOS of the exposed atoms (Ni-a, Ni-b, Ni-c, Ni-d, Ni-e, Ni-f, C-a, and C-b: see Fig. 4(a)) on the C-terminated surface are shown in Fig. S6, confirming that the surface states identified in Fig. 5(a) are unsaturated bonds between Ni and C, and that the states are split into two. Fig. S6 shows that the higher-energy surface states (ca. −3.9 eV) had a greater contribution from the three-coordinate carbon (C-b), whereas the lower-energy surface states (ca. −4.3 eV) exhibited a greater contribution from the five-coordinate carbon atom (C-a). The higher peak was attributed to the hybridisation between C-b and three neighbouring Ni atoms (Ni-a, Ni-c and Ni-e) (Fig. S6). Among the Ni atoms, Ni-c was the most abundant. Ni-c did not interact with the other exposed C atoms, in contrast to Ni-a and Ni-e. The different features in the PDOS are shown in Fig. S6. The projection to the magnetic angular momentum component confirms that the higher peak includes a p_z component from C-b (the z-direction is perpendicular to the surface) (Fig. S7(b)). The p_x and p_y components of C-b are present in the lower peak (Fig. S6(d) and S7(b)). In this lower-energy state, the contribution of the five-coordinate carbon, C-a, is larger (Fig. S6(d)) than that of C-b. Fig. S7(a) confirms that the p_z component of C-a is dominant in this peak. Hybridisation between C-a and the four neighbouring Ni atoms (Ni-a, Ni-b, Ni-e, and Ni-f) was detected in the PDOS, as shown in Fig. S6, whereas Ni-d, which did not bind to the exposed C atoms, exhibited a different DOS spectrum and no states in the energy range of the surface states. From these results, it can be concluded that Ni–C covalent interactions characterise the surface state of C-terminated Ni₃C(113).

To further prove the importance of Ni–C covalent interactions in the surface state, the state of the Ni-terminated surface was analysed. The PDOS of the surface-exposed atoms (Ni-g, Ni-h, Ni-i, Ni-j, Ni-k, Ni-l, and C-c; see Fig. 4(b)) on the Ni-terminated surface are shown in Fig. S8. The surface states generated at −4.3 eV are identified as hybridisation orbitals of Ni-g, Ni-h, Ni-k, and Ni-l with C-c, whereas the Ni-i and Ni-j states, which do not bond to the C-c atom, are not identified around −4.3 eV. Therefore, this surface state is accompanied by electron localisation on the Ni–C bonds. The projection onto the magnetic angular momentum component confirmed that the main contribution was from the p_z orbital (Fig. S7(c)). Thus, Ni–C covalent interactions are important, even on Ni-rich surfaces.

Owing to the importance of Ni–C covalent interactions in the surface states, the electronic and geometric structures can be explained by the rules and resonance structures presented in Section 3.1. First, surface-exposed carbon atoms, which will be characteristic features of carbides, are discussed. The resonance structure of three-coordinate C, C-b, was obtained by eliminating one boundary in Fig. 3(b) (Fig. 6, top panels). In this resonance structure, C-b has a formal dangling bond (an unpaired electron). To reduce the contribution of this configuration to the unpaired electrons, a double-bonding interaction is formed (Fig. 6, lower panel). Owing to the contributions of the resonance structures with Ni=C double bonds, the Ni–C-b bond length was significantly shorter than that in bulk Ni₃C (1.8838 Å) (Fig. 7(b)).

Fig. 6 Resonance structure of Ni₃C(113) with the three-coordinate C atom.

Fig. 7 Bond length between exposed C atoms in Ni₃C(113) and neighbouring Ni atoms: (a) Ni–C-a bonds, (b) Ni–C-b bonds, and (c) Ni–C-c bonds.

The resonance structures of the Ni atoms and five-coordinate carbons (C-a and C-c) are shown in Fig. 8. The topologies of C-a and C-c differ in terms of whether all the Ni atoms bound to the C atoms participate in other Ni–C bonds. When all the neighbouring Ni atoms of the C atoms are bonded to other C atoms, such as the C-a atom in the C-terminated surface, the resonance structure shown in Fig. 8(a) is obtained. The Ni atoms coordinated to C-a have a tetrahedral structure. In this arrangement, the four Ni atoms exposed on the surface are almost equivalent, whereas the Ni atoms in the axial position are not equivalent. Therefore, the contribution of the axial Ni (greyed structure in Fig. 8(a)) will be lower than that of the other resonance structures. A decrease in the resonance structure implies a decrease in the bond order. Indeed, Fig. 8(a) shows that the bond between the axial Ni and C is clearly long. When we assume that the contribution involving the four Ni atoms is dominant, the bond order is 3/6, and the formal charge of Ni is +6/6 (in fact, there is a small contribution from the axial position, so these values are upper limits). Both are larger than the bulk values (bond order: 1/3, formal charge: +2/3). The bond order and formal charge results derived from this resonance structure were consistent with the bond lengths shown in Fig. 7(a) and the charge analysis results in Table S3. The Ni–C bond lengths around C-a (except for axial Ni) are shorter than those in the bulk, and the surface-exposed Ni is more cationic.

Fig. 8 Resonance structures of Ni₃C(113) with five-coordinate C atoms (a) in C-terminated Ni₃C(113) and (b) in Ni-terminated Ni₃C(113).

Among the Ni atoms bonded to C-c on the Ni-terminated surface, two Ni atoms did not bond with the other C atoms. In this case, as shown in the resonance structure in Fig. 8(b), the contribution of the double bond between Ni and C is enhanced. Indeed, the Ni–C bond lengths around C-c are not symmetrical, and the two Ni–C-c bonds are shorter than those in bulk Ni₃C (Fig. 7(c)). Resonance structures in which the double bonds are cleaved are also possible; however, these structures involve two unpaired electrons. Hence, their contributions are little, and it can be assumed that no unpaired electrons were present in C-c (Fig. 8(b), top panels). Thus, for the five-coordinate C atoms, the electronic states are described by a resonance structure without formal dangling bonds.

The three-coordinate C atom on the C-terminated surface will be unstable because it has a resonance structure with states involving formal dangling bonds or a structure, in which the C–Ni–C bridging structure is broken. Therefore, it is assumed that a three-coordinate C atom translates into a five-coordinate carbon, for which the resonance structure can be described without the contribution of these unstable states. In our previous study,¹⁶ the structural optimisation of H adsorbed on Ni₃C(113) without symmetry restriction resulted in a significant surface deformation; in the optimised structure, the C-terminated surface no longer contained any three-coordinate carbons. Therefore, in this study, we performed structural optimisation using a surface structure without three-coordinate carbons as the initial structure. A reconstructed surface (Fig. 4(c)), in which all the C atoms exposed on the surface are five-coordinate C atoms, was converged. The estimated surface energy of the reconstructed surface was 0.220 eV Å⁻², which was 0.004 eV Å⁻² more stable than that of the C-terminated surface. Fig. S9 (SI) shows the DOS of the reconstructed surfaces. Similar to other surfaces, surface states with Ni–C correlations were identified. Therefore, the electronic state of the reconstructed surface follows the resonance structure of five-coordinate carbons, as shown in Fig. 8. Thus, using the resonance structure and the suggested rule based on electronic configuration, we predicted the surface stability of Ni₃C(113).

3.3. Hydrogen atom adsorption on Ni₃C(113) surfaces

In this section, we demonstrate that the theory based on the resonance structure along with the rule presented in Section 3.1 can explain the stability of H adsorption structures. For this purpose, we introduced the adsorption energy (E_ads), distortion energy (E_dis), and interaction energy (E_int).

H adsorption on C-terminated surfaces was reported in our previous study,¹⁶ but surface reconstruction was not considered. Therefore, the adsorption energy E_ads was recalculated based on the reconstructed surface. E_ads was calculated using eqn (1) to investigate whether molecular hydrogen is adsorbed by dissociation (i.e. dissociative adsorption is exothermic when E_ads < 0).

E_ads = E(ads/surf) − E(surf) − 0.5E(H₂) (1)

where E(H/surf) is the total energy of the H-adsorbed surface after structural optimisation. E(surf) is the total energy of the Ni₃C slab prior to H adsorption. The total energy of the reconstructed surface shown in Fig. 4(c) was used for the C-terminated surface, and the total energy of the slab model shown in Fig. 4(b) was used for the Ni-terminated surface. E(H₂) is the total energy of H₂. Since the results showed that Ni₃C(113) could be largely reconstructed, the surface distortion energy (E_dis) caused by H adsorption was calculated using eqn (2):

E_dis = E(surf)_fix − E(surf) (2)

where E(surf)_fix is the total energy of the slab models having the geometry after hydrogen adsorption (with the hydrogen atoms removed and the slab geometry fixed). E_int, defined by eqn (3), was used to investigate the strength of the interaction between hydrogen atoms and the surface.

E_int = E(H/surf) − E(surf)_fix − E(H) (3)

where E(H) is the total energy of the H atom.

The structures of the C-terminated surfaces with adsorbed H are denoted as 1-C, 2-C, 3-C, …, in decreasing order of stability and are discussed based on the results for the structures with H interacting with C atoms, which are shown in Fig. 9. When interacting with C, H was adsorbed at the top site of the C atom. The adsorption of C resulted in the most stable adsorption structure (1-C) and the least stable structure (10-C). In the 1-C structure, H interacts with a three-coordinate carbon atom (C-b), whereas in the 10-C structure, H interacts with a five-coordinate carbon atom (C-a). This indicates that the three-coordinate carbon adsorbed hydrogen more strongly than the five-coordinate carbon. This difference arises because (1) hydrogen adsorption on the five-coordinate carbon induces a large surface distortion (destabilisation), and (2) the interaction between the three-coordinate carbon and hydrogen is strong. Mechanism (1) is confirmed using the E_dis data presented in Fig. 9. Notably, among the Ni–C bonds with five-coordinate C atoms, the bond with Ni at the axial site is significantly elongated (Fig. 13(c), centre panel), resulting in a positive and large E_dis. Mechanism (2) was confirmed using the PDOS data (Fig. 10). As shown Fig. 10, the s state of H can be found in the energy region between −7.5 and −5.0 eV. This peak overlapped with that of the C p_z orbital. In the other regions, the DOS of H was small, and the partial charge distribution was small (Fig. S10 in the SI). The results of pCOHP analysis for the C–H bond shown in Fig. S11 indicate negative pCOHP below the Fermi energy and positive pCOHP above the Fermi energy for all structures. This is characteristic of typical covalent interactions. From these results, it is concluded that the C–H interactions are covalent bonds between the C p_z and H s orbitals.

Fig. 9 Geometry, adsorption energy (E_ads), distortion energy (E_dis), and interaction energy (E_int) of H adsorbed on the C atoms of C-terminated Ni₃C(113): (a) C-1, (b) C-5, and (c) C-10 structures.

Fig. 10 PDOS of H/C-terminated Ni₃C(113). (a) s component of H in 1-C, (b) p components of C in 1-C, (c) s component of H in 10-C, and (d) p components of C in 10-C. Positive and negative values represent majority and minority spins, respectively.

The PDOS (Fig. 10) and partial electron density distributions (Fig. S10) show peaks in the range of −7.5 to −5.0 eV, corresponding to σ bonding. Because the C–H interaction in H/C-terminated Ni₃C(113) is covalent, the strength of the interactions can be estimated from the shift in the PDOS of C p_z before (Fig. S7(a) and (b)) and after (Fig. 10(b) and (d)) H adsorption (i.e., a shift of 2.8 eV for the three-coordinate C atom (C-b) and 2.4 eV for the five-coordinate C atom (C-a)). Therefore, the three-coordinate C atom interacted more strongly with H, as indicated by the difference in E_int. This can also be confirmed from the integral values of pCOHP, specifically 1-C structure: −3.095 eV, 5-C structure: −3.106 eV, and 10-C structure: −3.041 eV. The order of pCOHP integral values corresponds to the order of E_int values, indicating that the interaction between three-coordinate C and H is stronger than that between five-coordinate C and H, although the difference is not large. Comparing the results of E_int and E_dis (Fig. 9), the energy difference between the 1-C and 10-C adsorption structures is larger for E_dis; hence, it is assumed that mechanism (1), that is, surface distortion, is the main reason for the difference in the E_ads values of C-a and C-b. This is supported by the pCOHP results shown in Fig. S12. The shapes of the pCOHP plots for C–Ni pairs exposed on the surface are very similar, and their peak intensities correspond well with the order of E_ads and E_dis. The specific integral values are −1.433 eV (1-C), −1.395 eV (5-C), and −1.100 eV (10-C). Furthermore, the pCOHP between the axial Ni and C in the 10-C structure is exceptionally small compared to the others. Its integral value is −0.1939 eV: it is one order of magnitude smaller than the others. This pCOHP results indicate that the interaction with the axial Ni is significantly weakened by H adsorption. Similar to the 1-C structure, H was adsorbed on the three-coordinate C atom in the 5-C structure. However, the surface distortion was too large, and the stabilisation due to adsorption was not significant (Fig. 9). The Ni atoms around the three-coordinate C atom in the 1-C structure are arranged in an almost equilateral-triangular structure, whereas those in the 5-C adsorption structure are arranged in a significantly distorted arrangement (Fig. S13 in the SI). The same conferment was obtained from pCOHP (Fig. S12). Based on the integral values, it can be concluded that the difference between the 1-C and 5-C structures is largely attributable to the changes in the covalent character of the Ni–C bond.

For the H/Ni-terminated Ni₃C(113) system, the optimised structures were denoted as 1-Ni, 2-Ni, and 3-Ni, starting from the most stable adsorbed structure. The Ni-terminated (113) surface had only one surface-exposed carbon (C-c) per calculated unit cell and only one structure in which the adsorbed H interacted with C. The E_ads, E_dis, E_int, and PDOS of C and H and the partial electron density distribution of the C–H binding region are summarised in Fig. 11. The structure shown in Fig. 11 was the least stable for H adsorption on the Ni-terminated surface (6-Ni). In this adsorption structure, H interacts with the five-coordinate C (C-c), and the C–H bond is covalent, as confirmed by the PDOS and partial charge density. With the formation of a C–H covalent bond, the bond between C and axial Ni was elongated, and E_dis was positive and large. Because of this large surface destabilisation, E_ads was positive (the dissociative adsorption of molecular hydrogen is thermally unstable). These results are similar to those obtained for H–C-a interactions in the 10-C structure. As indicated by the E_dis values, H adsorption on C-c on the Ni-terminated surface (Fig. 11) had a lower energy than that on C-a (Fig. 9(c)). This is because the elongation of the distance between the axial Ni and the five-coordinate carbon atoms is larger for the C-terminated surface.

Fig. 11 Calculated results for H adsorption on the C on the top-site of Ni-terminated Ni₃C(113) (6-Ni structure). (a) Geometry, adsorption energy, surface distortion energy, and H–Ni₃C interaction energy. (b) PDOS of the H atom. (c) PDOS of the C atom interacting with H. (d) Partial charge density distribution in the energy range between −7.5 and −5.0 eV (the threshold is 0.01 a.u.).

The difference in the H adsorption onto the C atoms can also be explained by the resonance structures of Ni₃C. As shown in Fig. 6, the three-coordinate C atoms contain unpaired electrons as dangling bonds. In other words, the valence state of the C atom is unity. The valence state of H is also unity. Hence, the three-coordinate C and H atoms formed a strong covalent bond, and the C–H interaction resulted in H₂ dissociation. The five-coordinate C atoms had no formal unpaired electrons (full-fill electron configuration), as shown in Fig. 8. Therefore, one of the Ni–C bonds must be broken to form a bond with H, resulting in a dangling bond. The axial Ni, which has the longest Ni–C bond, dissociates and leaves C with a valence of unity, after which the C atom bonds with hydrogen. Because of this Ni–C bond cleavage, E_dis becomes too large for the adsorption of H on the five-coordinate carbon, making it difficult to stabilise the adsorbed state sufficiently for H–H bond cleavage. Specifically, the resonance structures of H-adsorbing three- and five-coordinate C atoms were drawn, as shown in Fig. 12. These resonance structures were written according to the rules presented in Section 3.1. As shown by the resonance structures, H adsorption on the three-coordinated carbon did not weaken the Ni–C bond (the bond order was 1/3), whereas H adsorption on the five-coordinate carbon reduced the bond order of the Ni–C bond. The pCOHP results shown in Fig. S11 and S12 are also consistent with the resonance structure in Fig. 12. The resonance structure in Fig. 12(b) indicates that H adsorption on the five-coordinate carbon caused a large surface distortion and weakened the C–H interaction. In addition, we can explain the difference between the C-terminated and Ni-terminated surfaces. The five-coordinate C on the Ni-terminated surface was adjacent to the Ni atoms that were not bonded to any other C atoms, and a double bond was formally introduced between C and Ni (Fig. 8(b)). The five-coordinate carbon on the Ni-terminated surface can also generate unpaired electrons by converting this double bond to a single bond instead of severing the bond with the Ni immediately below. Therefore, as shown in Fig. 11, the change in the bond length between the five-coordinate C and axial C in the H/Ni-terminated Ni₃C was smaller than that in the H/C-terminated Ni₃C, and the Ni–C distance with the formal double bond increased when H was adsorbed on the Ni-terminated Ni₃C. Thus, using the rule and resonance structure, we can predict the stability of H adsorption onto C atoms and structural changes without the specific results of adsorption energy, density of states, electron distribution, or geometry optimisation of H/Ni₃C systems.

Fig. 12 Resonance structures of H adsorption onto (a) 3-coordinate and (b) 5-coordinate C atoms.

Fig. S14 and S15 (SI) summarise the results of the converged structures, in which H interacts with Ni atoms. The structures, in which H was adsorbed on the hollow sites surrounded by three Ni atoms, were stable, and the Ni–H interaction was stronger on the Ni-terminated surface. Adsorption onto the Ni–Ni bridge sites also converged; however, these interactions were less stable than those of adsorption onto the hollow sites. Fig. S16 shows the PDOS and partial electron-density distributions of the 2-C and 1-Ni structures. Although Ni₃C has a Ni d band (>−4 eV), the states of the H atoms are located in the region between −7.5 and −4.0 eV, where the states originating from C atoms have a larger distribution than the states originating from Ni atoms (see Fig. S16). Therefore, even after the adsorption of H on Ni, the contribution of C cannot be neglected in Ni₃C. As a clear result of the contribution of C, we show that E_ads for H adsorption onto the Ni hollow sites in Ni₃C(113) follows a trend dictated by the number of Ni–C bonds of Ni interacting with H (N_C–Ni). Fig. 13 shows a clear correlation between E_ads and N_C–Ni, where E_ads became more negative as N_C–Ni decreased. The correlation between the E_ads of N_C–Ni can be explained by the formal Ni–C double bond in the resonance structure. Ni with two bonds to C in the topology shown in Fig. 3 has a formal valence of zero, whereas Ni with one bond can assume a valence of unity, as shown in Fig. 8. Hollow sites with less N_C–Ni can adsorb H more strongly. The results shown in Fig. 13 and Fig. S16 indicate that Ni–C covalent interactions are necessary for understanding Ni–H interactions and the importance of the resonance structure.

Fig. 13 Correlation between adsorption energy (E_ads) and summation of the number of neighbouring C atoms of adsorbed Ni (N_C–N).

The Ni₃C–H interaction can be rationalised through resonance structures, suggesting that the charge on the adsorbed H atom correlates with electronegativity differences. Table 2 lists the Bader charges of H and interacting atoms in various adsorption configurations. When the H atom interacts with Ni atoms in Ni₃C, it carries a negative charge, which is consistent with the interactions between Ni metal and H atoms.¹⁶ In contrast, the H atoms interacting with C atoms are almost neutral or slightly cationic. This result is consistent with the fact that H has a slightly lower electronegativity than C and that in molecules composed only of C and H, such as CH₄, the electron pairs are slightly skewed towards C. The formation of polar hydrogen species on the catalyst surface is beneficial for hydrogenation reactions.¹⁷ A surface structure that generates polar hydrogen species, that is, a surface with exposed Ni and C species, as determined by covalent bonding, is crucial for understanding the catalytic activity of Ni₃C.

Table 2 Atomic charges of H in H/Ni₃C(113) systems

Adsorption structure	Adsorbing atom	Charge of adsorbed H [e^−]
1-C	C	+0.05
2-C	Ni	−0.28
3-C	Ni	−0.27
4-C	Ni	−0.30
5-C	C	+0.10
6-C	Ni	−0.25
7-C	Ni	−0.32
8-C	Ni	−0.32
9-C	Ni	−0.24
10-C	C	+0.08
1-Ni	Ni	−0.36
2-Ni	Ni	−0.33
3-Ni	Ni	−0.31
4-Ni	Ni	−0.29
5-Ni	Ni	−0.29
6-Ni	C	+0.06

---

3.4. Results and discussion on other metal carbides

The objective of this study is to provide a theoretical framework beyond the d-band centre theory for understanding the unique hydrogen activation property of Ni₃C. As previously discussed, in Ni₃C, hybridised carbon orbitals have not been confirmed; instead, bands attributable to ionic bonding, covalent bonding, and metallic bonding have been confirmed. The surface levels of Ni₃C are significantly influenced by the covalent bands, and their electronic states and interactions with hydrogen can be understood through resonance structures. This section briefly discusses whether the insights from Ni₃C can be extended to understanding the catalytic activities of other metal carbides. The calculated results of Co₂C, Fe₂C, Fe₃C, CrC, Mo₂C, and WC are summarised in the SI (Fig. S17). The calculation conditions were identical to those for Ni₃C; electronic state and structural optimisation was performed using the structures obtained from the Materials Project as the initial structures.

The PDOS of Co₂C is similar to that of Ni₃C. This can be understood from the electronic configuration of Ni₃C shown in Fig. 3(a). Since Co²⁺ has one fewer electron than Ni²⁺, it adopts a configuration where one d electron is removed from Ni²⁺: namely, a d⁷ state with three unpaired valence electrons. In the Co₂C structure shown in Fig. 14(a), Co has three bonds to bridge C atoms. That is, while Ni bridges two C atoms, Co bridges three C atoms. Consequently, Co₂C can be depicted with a similar resonance structure to Ni₃C, and its electronic state is analogous to that of Ni₃C. Fe₂C involves a further reduction of one or more d electrons from the metal species compared to Co₂C. Since Fe²⁺ is a d⁶ electron system and Fe³⁺ is a d⁵ electron system, spin polarisation is expected to occur in the d electrons not utilised for covalent bonding with carbon. Nevertheless, no significant hybridisation between the s and p orbitals of carbon is detected, which is consistent with Ni₃C and Co₂C. In addition, the covalent interaction and metallic Fe band can be found in both iron carbides. Although Cr, Mo, and W are in the same group, the DOS of carbides shows marked differences. CrC has a similar electronic structure to Ni₃C, Co₂C, Fe₂C, and Fe₃C: CrC has a Cr metallic band on its Fermi level and a metal–carbon correlation band below the metallic band. Compared to Ni₃C and Co₂C, the boundary between the ionic bonding band and the covalent bonding band is not clearly defined in the Mo₂C. The PDOS of WC confirms that its band gap is wider than that of Ni₃C, and the absolute value of the formal charge of carbon is clearly greater than 2. Therefore, interactions with a larger contribution from ionic bonding than covalent bonding may arise in the WC. However, the degree of electron localisation (ionisation) depends on the on-site Coulomb correction. These results indicate that the balance between ionic, covalent, and metallic bonding will depend on the metal species.

Fig. 14 Comparison of results calculated using DFT and predicted by PFP. (a) E_ads for H/C-terminated Ni₃C(113), (b) E_dis for H/C-terminated Ni₃C(113), (c) E_int for H/C-terminated Ni₃C(113), (d) number of valence electrons of adsorbed H for H/C-terminated Ni₃C(113) estimated using the Bader algorithm, (e) E_ads for H/Ni-terminated Ni₃C(113), (f) E_dis for H/Ni-terminated Ni₃C(113), (g) E_int for H/Ni-terminated Ni₃C(113), and (h) number of valence electrons of adsorbed H on H/Ni-terminated Ni₃C estimated using the Bader algorithm (113). Regression results are summarised in Table S4 of the SI.

Thus, while there are some differences in the electronic structures of the calculated metal carbides, the following characters are common: no significant hybridisation is found between the s and p orbitals of carbon; orbital correlation is found between carbon and metal species (covalent character); the band with metal alone is found (metallic character); and charge transfer occurs between the metal and carbon (ionic character). The differences are just in the positions of the bands originating from these three characters and the magnitude of their contributions. Determining this correct balance will require a detailed discussion of the electron exchange term. As shown in Fig. S2, the band gap between the covalent bonding band and the metallic bonding band disappears as the on-site Coulomb parameter is increased (taking into account the repulsion between d electrons and making magnetisation easier). Hence, for establishment of general rules for the catalytic activity of metal carbides, it is necessary to correctly understand the influence of the electron exchange term compared with experiments. However, based on the common features depicted above, it can be concluded that the contribution of the covalent bond band, as demonstrated by the resonance structure in this study, will prove useful in discussions of the catalytic activity of other metal carbides, albeit to varying degrees.

Towards systematising the catalytic activity of metal carbides, numerous other metal carbides and their surface models should be investigated theoretically. Hence, the use of computational methods faster than DFT calculations will be needed to establish general rules. In the next section, we discuss whether the MLIP calculation is effective for theoretical investigations, which is 100–1000 times faster than DFT.

3.5. Comparison of DFT and MLIP

The significance of resonance structures in elucidating the interaction between Ni₃C and H, particularly in explaining the valence charge distribution and adsorption energies of H on Ni₃C surfaces, has been demonstrated above. The C in Ni₃C does not exhibit sp³ hybridisation; instead, C–H bonding primarily occurs via the p_z orbital. Additionally, surface reconstructions, in which three-coordinate C atoms moved to a five-coordinate site, were observed for the C-terminated surfaces. H adsorption further induces pronounced surface distortions and modifies the associated resonance structure. In this section, we investigate the feasibility of employing a universal MLIP for surface reactions, specifically the PFP model,⁴⁹ to predict Ni₃C–H interactions based on the framework of resonance structure theory.

Fig. 14 presents a comparison of the DFT and PFP data; the results for H/C-terminated Ni₃C(113) are shown in Fig. 14(a)–(d), and those for H/Ni-terminated Ni₃C(113) are shown in Fig. 14(e)–(h); comparisons of the E_ads, E_dis, E_int, and the valence electrons of H estimated using the Bader algorithm are shown in Fig. 14(a), (e), (b), (f), (c), (g), (d) and (h), respectively. The DFT–MLIP correlations of E_ads for the C-terminated surface (Fig. 14(a)) revealed a strong agreement overall, with one notable outlier (indicated by red dots). This outlier corresponds to the structure where the H atom is adsorbed onto five-coordinate carbon atoms. As shown in Fig. 14(b) and (c), the errors in E_dis were particularly pronounced. The correlation between the E_ads data for the Ni-terminated surface obtained using DFT and MLIP was generally good (Fig. 14(e)); however, for H adsorption onto the five-coordinate C atoms, the MLIP predicted the wrong sign of E_ads (blue dot). Moreover, E_int for H adsorption on five-coordinate C obtained using PFP was clearly an outlier (blue dots in Fig. 14(g)). A previous study has shown that PFP is not well-suited for capturing the resonance destabilisation caused by molecular adsorption onto solids.⁷⁴ Bond lengths are good descriptors for estimating changes in resonance structures; the predicted structures with C–H bonds by PFP are summarised in Fig. 15. Remarkably, the bond distances in the structures without hydrogen adsorption (left panels of Fig. 15) are in good agreement with the calculated results by DFT (Fig. 7). The hydrogen adsorption structure onto the three-coordinate carbon is correctly predicted (Fig. 15(a)), whereas the adsorption structure onto the five-coordinate carbon of the C-rich surface (red dot in Fig. 14, right panel of Fig. 15(b)) is clearly incorrect. The bonding distance with the axial Ni is so long, and the quadrilateral formed by the other four Ni atoms is rhombic (a resonance structure where the rhombic shape is dominant cannot be written). To accurately estimate the surface distortion produced by the C-terminated surface, it is essential to account for changes in the resonance of the Ni–C bond shown in Fig. 12, particularly the significant weakening of one Ni–C bond. It is likely that PFP underestimates E_dis because of the underestimation of the decrease in stabilisation caused by the resonance of the Ni–C bond. For a system with a Ni-terminated surface, there was no significant difference in the structure before and after hydrogen adsorption structures estimated by PFP and DFT (Fig. 15(c)). However, understanding the C–H bond requires an understanding of the destabilisation caused by the decrease in the contribution of resonance structures with Ni–C double bonds. This resonance destabilisation would not be learned correctly—the reason for significant overestimation of E_int the PFP model, as shown in Fig. 14(g).

Fig. 15 Bond lengths around surface-exposed carbon atoms in the optimised structures with (left panels) and without (right panels) H adsorption onto the on-top sites of carbons by PFP: (a) three-coordinate carbon in C-rich Ni₃C(113), (b) five-coordinate carbon in C-rich Ni₃C(113), and (c) five-coordinate carbon in Ni-rich Ni₃C(113).

Although some challenges remain in accurately estimating the resonance destabilisation associated with five-coordinate C, PFP correctly predicts the most stable H adsorption configuration. Furthermore, PFP can predict the charge of the adsorbed H atoms, and the estimated values are qualitatively correct (Fig. 14(d) and (h)): H on Ni atoms is anionic, and H on C atoms is neutral or slightly cationic. These results support the conclusion that PFP is sufficiently accurate for a qualitative discussion of hydrogen activation by Ni₃C, whereas careful analysis of the mechanism involving five-coordinate C atoms is necessary.

4. Conclusion

We propose that resonance structure theory is effective for understanding the interaction of H with Ni₃C, which has high catalytic activity in hydrogen evolution and hydrogenation reactions. The resonance structure was drawn according to the following rules:

(1) Formal charge and valence of C are −2 and 2, respectively;

(2) Bonding of Ni to C²⁻ results in a formal charge of +2; and

(3) The formal valence of Ni²⁺ is 2.

The nature of the surface states and stability of the adsorbed H can be inferred from the resonance structure based on the rules discussed above. Specifically, the following conclusions have been drawn: (A) three-coordinate carbons possess formal, unpaired electrons (dangling bonds); (B) surface reconstruction occurs on C-terminated surfaces to reduce the contribution of these dangling bonds to the three-coordinate carbons; (C) H adsorption on three-coordinate carbon yields the most stable configuration due to the presence of dangling bonds; (D) hydrogen adsorption on five-coordinate carbon leads to unstable configurations primarily due to significant weakening of the Ni–C bond; (E) the adsorption energies of H onto the hollow Ni sites are correlated with the coordination number of Ni to C; (F) the formal charge of adsorbed hydrogen is governed by its electronegativity relative to the host atom; and (G) hydrogen on Ni₃C undergoes heterolytic dissociation, which is beneficial for hydrogenation reactions.

To derive these insights, it is essential to correctly estimate changes in the resonance structure. However, even PFP, the most accurate universal potential for surface reactions available at the time of this study, falls short in capturing these resonance effects adequately. Therefore, developing algorithms that can explicitly learn and represent resonance structures remains a key challenge for MLIP approaches. Despite these limitations, PFP successfully predicted the most stable adsorption configurations and heterogeneous charges of the adsorbed hydrogen atoms (H⁻ on Ni and H⁰/H⁺ on C). This indicates that PFP is a useful tool for screening hydrogen adsorption structures and predicting their associated charges. The rule for the resonance structure of Ni₃C presented in this study will be important for further understanding the properties of metal carbides. Further validation of the application range of the rules and theories presented in this study is an important direction for future research. A theoretical investigation in this regard is currently underway with acceleration by MLIP.

Author contributions

K. T., S. Y., and T. M. conceived the study. K. T. designed the theoretical calculations. K. T. and K. M. performed the calculations. K. T. and K. M. analysed the results with inputs from S. Y., R. K., Y. K., and T. M. R. K. and Y. K. validated the calculated results. K. T. wrote the original draft and all authors commented on the manuscript.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d5cp02869h. The contents of SI are as following: (1) dependence on k-points sampling, (2) results by DFT+U method, (3) comparison of DFT and CHGNet, (4) PDOS of Ni and C in bulk Ni₃C, (5) results for dual-sided slab models, (6) PDOS results for the C-terminated Ni₃C(113), (7) PDOS results for the C atoms in Ni₃C(113), (8) PDOS results for Ni-terminated Ni₃C(113), (9) DOS and PDOS results for the reconstructed Ni₃C(113), (10) partial electron density distribution of 1-C and 10-C, (11) pCOHP results, (12) comparison of 1-C and 5-C, (13) optimised structures of the H/C-terminated surfaces, (14) optimised structures of the H/Ni-terminated surfaces, (15) electronic states of H in the hollow sites of Ni in the Ni₃C, (16) PDOS of other several metal carbides, (17) comparison of population analysis methods, (18) lattice parameters of Ni₃C, (19) atomic charges of exposed atoms in Ni₃C(113), (20) regression results for DFT–MLIP correlation, (21) estimation scheme of surface energy, and (22) specific coordinates of optimised structures.

Acknowledgements

The computations in this study were performed using the facilities at the Supercomputer Centre, Institute for Solid State Physics, University of Tokyo (project code: 2025-Ca-0013), and the Research Institute for Information Technology, Kyushu University. This study was conducted under the auspices of the Japan Society for the Promotion of Science (JSPS KAKENHI; grant numbers JP22H02050, JP22H04974, JP24H00459, JP24K01255, and JP25K08589) and the Japan Science and Technology Agency (JST) PRESTO program (grant number JPMJPR24MB).

Notes and references

1. B. Hammer and J. K. Nørskov, Nature, 1995, 376, 238–240,  DOI:10.1038/376238a0.
2. V. Stamenkovic, B. S. Mun, K. J. J. Mayrhofer, P. N. Ross, N. M. Markovic, J. Rossmeisl, J. Greeley and J. K. Nørskov, Angew. Chem., Int. Ed., 2006, 45, 2897–2901,  DOI:10.1002/anie.200504386.
3. A. Kulkarni, S. Siahrostami, A. Patel and J. K. Nørskov, Chem. Rev., 2018, 118, 2302–2312,  DOI:10.1021/acs.chemrev.7b00488.
4. W. S. Williams, Prog. Solid State Chem., 1971, 6, 57–118,  DOI:10.1016/0079-6786(71)90028-8.
5. J. S. Lee, in Encyclopedia of Catalysis, ed. I. Horváth, John Wiley & Sons, Chichester, 2010 DOI:10.1002/0471227617.eoc136.pub2.
6. L. Borchardt, C. Hoffmann, M. Oschatz, L. Mammitzsch, U. Petasch, M. Herrmann and S. Kaskel, Chem. Soc. Rev., 2012, 41, 5053–5067,  10.1039/C2CS15324F.
7. A. M. Alexander and J. S. J. Hargreaves, Chem. Soc. Rev., 2010, 39, 4388–4401,  10.1039/B916787K.
8. H. Wang and L. Gao, Curr. Opin. Electrochem., 2018, 7, 7–14,  DOI:10.1016/j.coelec.2017.10.010.
9. W. Yang, S. Rehman, X. Chu, Y. Hou and S. Gao, ChemNanoMat, 2015, 1, 376–398,  DOI:10.1002/cnma.201500073.
10. X. Fan, Z. Peng, R. Ye, H. Zhou and X. Guo, ACS Nano, 2015, 9, 7407–7418,  DOI:10.1021/acsnano.5b02420.
11. H. Wang, Y. Cao, G. Zou, Q. Yi, J. Guo and L. Gao, ACS Appl. Mater. Interfaces, 2017, 9, 60–64,  DOI:10.1021/acsami.6b14393.
12. H. Fan, H. Yu, Y. Zhang, Y. Zheng, Y. Luo, Z. Dai, B. Li, Y. Zong and Q. Yan, Angew. Chem., Int. Ed., 2017, 56, 12566–12570,  DOI:10.1002/anie.201706610.
13. J. G. Roblero, F. Pola-Albores, M. A. Valenzuela, E. Rojas-García, E. Ríos-Valdovinos and G. Valverde-Aguilar, Int. J. Hydrogen Energy, 2019, 44, 10473–10483,  DOI:10.1016/j.ijhydene.2019.02.119.
14. R. F. André, L. Meyniel and S. Carenco, Catal. Sci. Technol., 2022, 12, 4572–4583,  10.1039/D2CY00894G.
15. K. He, J. Xie, Z.-Q. Liu, N. Li, X. Chen, J. Hu and X. Li, J. Mater. Chem. A, 2018, 6, 13110–13122,  10.1039/C8TA03048K.
16. S. Yamaguchi, D. Kiyohira, K. Tada, T. Kawakami, A. Miura, T. Mitsudome and T. Mizugaki, Chem. – Eur. J., 2024, 30, e202303573,  DOI:10.1002/chem.202303573.
17. T. Kawakami, S. Yamaguchi, S. Suganuma, K. Nakajima, T. Mitsudome and T. Mizugaki, ACS Sustainable Chem. Eng., 2025, 13, 7994–8002,  DOI:10.1021/acssuschemeng.5c01806.
18. R. Shen, K. He, A. Zhang, N. Li, Y. H. Ng, P. Zhang, J. Hu and X. Li, Appl. Catal., B, 2021, 291, 120104,  DOI:10.1016/j.apcatb.2021.120104.
19. Y. Hirayama, A. Miura, M. Hirayama, H. Nakamura, K. Fujita, H. Kageyama, S. Yamaguchi, T. Mizugaki and T. Mitsudome, Small, 2025, 21, 2412217,  DOI:10.1002/smll.202412217.
20. S. Yamaguchi, M. Sarmiento-Diaz, T. Kawakami, K. Tada, T. Mitsudome and T. Mizugaki, ACS Appl. Nano Mater., 2025, 8, 13886–13893,  DOI:10.1021/acsanm.5c02582.
21. Y. Liu, T. G. Kelly, J. G. Chen and W. E. Mustain, ACS Catal., 2013, 3, 1184–1194,  DOI:10.1021/cs4001249.
22. J. R. Kitchin, J. K. Nørskov, M. A. Barteau and J. G. Chen, Catal. Today, 2005, 105, 66–73,  DOI:10.1016/j.cattod.2005.04.008.
23. D. Tian, S. R. Denny, K. Li, H. Wang, S. Kattel and J. G. Chen, Chem. Soc. Rev., 2021, 50, 12338–12376,  10.1039/D1CS00590A.
24. S. Posada-Pérez, F. Viñes, J. A. Rodriguez and F. Illas, Top. Catal., 2015, 58, 159–173,  DOI:10.1007/s11244-014-0355-8.
25. M. G. Quesne, A. Roldan, N. H. de Leeuw, C. R. A. Richard and A. Catlowa, Phys. Chem. Chem. Phys., 2019, 21, 10750–10760,  10.1039/C9CP00924H.
26. V. Ellingsson, A. Iqbal, E. Skúlason and Y. Abghoui, ChemSusChem, 2023, 16, e202300947,  DOI:10.1002/cssc.202300947.
27. L. Wang, B. Wang, M. Fan, L. Ling and R. Zhang, Fuel, 2023, 336, 127131,  DOI:10.1016/j.fuel.2022.127131.
28. X. Tang, R. Salehin, G. B. Thompson and C. R. Weinberger, Phys. Rev. Mater., 2021, 5, 103603,  DOI:10.1103/PhysRevMaterials.5.103603.
29. M. López, F. Viñes, M. Nolan and F. Illas, J. Phys. Chem. C, 2020, 124, 15969–15976,  DOI:10.1021/acs.jpcc.0c03893.
30. X. Shen, D. Luo, J. Ni, J. Wu, C. Ma, H. Suo, X. Liu, Z. Lv, X. Wen, Y. Li, T. Zhang and Y. Yang, Appl. Surf. Sci., 2021, 570, 151210,  DOI:10.1016/j.apsusc.2021.151210.
31. C. Kunkel, F. Viñes and F. Illas, J. Phys. Chem. C, 2019, 123, 3664–3671,  DOI:10.1021/acs.jpcc.8b11942.
32. R. Zhang, L. Kang, H. Liu, L. He and B. Wang, Comput. Mater. Sci., 2018, 145, 263–279,  DOI:10.1016/j.commatsci.2018.01.013.
33. E. Osei-Agyemang, J.-F. Paul, R. Lucas, S. Foucaud and S. Cristol, J. Phys. Chem. C, 2016, 120, 8759–8771,  DOI:10.1021/acs.jpcc.6b01460.
34. X.-W. Liu, C.-F. Huo, Y.-W. Li, J. Wang and H. Jiao, Surf. Sci., 2012, 606, 733–739,  DOI:10.1016/j.susc.2011.12.018.
35. P.-P. Chen, J.-X. Liu and W.-X. Li, ACS Catal., 2019, 9, 8093–8103,  DOI:10.1021/acscatal.9b00649.
36. L. Lin, X. Yang, P. Shi, L. Yan, K. Xie, C. Deng and Z. Chen, Surf. Interfaces, 2023, 40, 103100,  DOI:10.1016/j.surfin.2023.103100.
37. N. Kosar, S. Munsif, K. Ayub and T. Mahmood, Int. J. Hydrogen Energy, 2021, 46, 9163–9173,  DOI:10.1016/j.ijhydene.2021.01.011.
38. R. Gubo, P. Ren, X. Yu, T. Zhang, X. Wen, Y. Yang, Y. W. Li, J. W. Niemantsverdriet and C. J. Weststrate, Appl. Surf. Sci., 2023, 626, 157245,  DOI:10.1016/j.apsusc.2023.157245.
39. F. Silveri, M. G. Quesne, A. Roldan, N. H. de Leeuw and C. R. A. Catlow, Phys. Chem. Chem. Phys., 2019, 21, 5335–5343,  10.1039/C8CP05975F.
40. B. Ding, W.-J. Ong, J. Jiang, X. Chen and N. Li, Appl. Surf. Sci., 2020, 15, 143987,  DOI:10.1002/chem.201804686.
41. B. Huang, N. Zhou, X. Chen, W.-J. Ong and N. Li, Chem. – Eur. J., 2018, 10, 18479–18486,  DOI:10.1002/chem.201804686.
42. Z. Zeng, X. Chen, K. Weng, Y. Wu, P. Zhang, J. Jiang and N. Li, npj Comput. Mater., 2021, 7, 80,  DOI:10.1038/s41524-021-00550-4.
43. J. Jiang, Y. Zou, A. Arramel, F. Li, J. Wang, J. Zou and N. Li, J. Mater. Chem. A, 2021, 9, 24195–24214,  10.1039/D1TA07332J.
44. Y. Ma, J. Liu, M. Chu, J. Yue, Y. Cui and G. Xu, Catal. Lett., 2020, 150, 1418–1426,  DOI:10.1007/s10562-019-03033-w.
45. S. Smed, T. Salmi, L. P. Lindfors and O. Krause, Appl. Catal., A, 1996, 144, 177–194,  DOI:10.1016/0926-860X(96)00103-2.
46. F. Hu, J. Peng, W. Xie and N. Li, RSC Adv., 2021, 12, 869–873,  10.1039/D1RA07448B.
47. C. Chen and S. P. Ong, Nat Comput. Sci., 2022, 2, 718–728,  DOI:10.1038/s43588-022-00349-3.
48. B. Deng, P. Zhong, K. Jun, J. Riebesell, K. Han, C. J. Bartel and G. Ceder, Nat. Mach. Intell., 2023, 5, 1031–1041,  DOI:10.1038/s42256-023-00716-3.
49. S. Takamoto, C. Shinagawa, D. Motoki, K. Nakago, W. Li, I. Kurata, T. Watanabe, Y. Yayama, H. Iriguchi, Y. Asano, T. Onodera, T. Ishii, T. Kudo, H. Ono, R. Sawada, R. Ishitani, M. Ong, T. Yamaguchi, T. Kataoka, A. Hayashi, N. Charoenphakdee and T. Ibuka, Nat. Commun., 2022, 13, 2991,  DOI:10.1038/s41467-022-30687-9.
50. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868,  DOI:10.1103/PhysRevLett.77.3865.
51. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979,  DOI:10.1103/PhysRevB.50.17953.
52. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775,  DOI:10.1103/PhysRevB.59.1758.
53. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558–561,  DOI:10.1103/PhysRevB.47.558.
54. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 14251–14269,  DOI:10.1103/PhysRevB.49.14251.
55. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186,  DOI:10.1103/PhysRevB.54.11169.
56. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50,  DOI:10.1016/0927-0256(96)00008-0.
57. H. J. Monkhorst and J. D. Pack, Phys. Rev. B: Condens. Matter Mater. Phys., 1976, 13, 5188–5192,  DOI:10.1103/PhysRevB.13.5188.
58. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104,  DOI:10.1063/1.3382344.
59. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465,  DOI:10.1002/jcc.21759.
60. V. I. Anisimov, I. V. Solovyev, M. A. Korotin, M. T. Czyzyk and G. A. Sawatzky, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 48, 16929,  DOI:10.1103/PhysRevB.48.16929.
61. G. Henkelman, A. Arnaldsson and H. Jónsson, Comput. Mater. Sci., 2006, 36, 354–360,  DOI:10.1016/j.commatsci.2005.04.010.
62. E. Sanville, S. D. Kenny, R. Smith and G. Henkelman, J. Comput. Chem., 2007, 28, 899–908,  DOI:10.1002/jcc.20575.
63. W. Tang, E. Sanville and G. Henkelman, J. Phys.: Condens. Matter, 2009, 21, 084204,  DOI:10.1088/0953-8984/21/8/084204.
64. M. Yu and D. R. Trinkle, J. Chem. Phys., 2011, 134, 064111,  DOI:10.1063/1.3553716.
65. V. L. Deringer, A. L. Tchougréeff and R. Dronskowski, J. Phys. Chem. A, 2011, 115, 5461–5466,  DOI:10.1021/jp202489s.
66. S. Maintz, V. L. Deringer, A. L. Tchougréeff and R. Dronskowski, J. Comput. Chem., 2013, 34, 2557–2567,  DOI:10.1002/jcc.23424.
67. R. Nelson, C. Ertural, J. George, V. L. Deringer, G. Hautier and R. Dronskowski, J. Comput. Chem., 2020, 41, 1931–1940,  DOI:10.1002/jcc.26353.
68. C. Ertural, S. Steinberg and R. Dronskowski, RSC Adv., 2019, 9, 29821–29830,  10.1039/C9RA05190B.
69. K. Momma, F. Izumi and VESTA, J. Appl. Crystallogr., 2011, 44, 1272–1276,  DOI:10.1107/S0021889811038970.
70. S. Nagakura, J. Phys. Soc. Jpn., 1958, 13, 1005–1014,  DOI:10.1143/JPSJ.13.1005.
71. J. Y. Hwang, A. R. P. Singh, M. Chaudhari, J. Tiley, Y. Zhu, J. Du and R. Banerjee, J. Phys. Chem. C, 2010, 114, 10424–10429,  DOI:10.1021/jp102571g.
72. C. M. Fang, M. H. F. Sluiter, M. A. van Huis and H. W. Zandbergen, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 86, 134114,  DOI:10.1103/PhysRevB.86.134114.
73. E. C. Stoner, Proc. R. Soc. London, Ser. A, 1938, 165, 372–414,  DOI:10.1098/rspa.1938.0066.
74. K. Tada, Chem. Lett., 2024, 53, upae173,  DOI:10.1093/chemle/upae173.

---