Electrochemical hydrogen evolution on Ni₂P: insights from constant-potential DFT and microkinetic modelling

Abstract

Hydrogen evolution on nickel phosphides (Ni₂P) is strongly influenced by surface structure and applied potential. Density functional simulations in combination with the constant-potential approach and a hybrid solvation model show that on Ni₂P(001), pristine Ni₃P₂ terminations favour Volmer–Volmer–Tafel energy landscapes, while phosphorus-reconstructed surfaces favour Volmer–Heyrovsky pathways. However, steady-state microkinetic modelling reveals that kinetic competition under bias shifts the predominant current response toward Volmer–Heyrovsky cycles on both terminations. The reconstructed surface achieves substantially lower overpotentials (120–180 mV at 10–100 mA cm⁻²), in close agreement with experimental measurements, while the pristine surface remains less active (330–370 mV). The rate is limited by Volmer on the pristine surface and Heyrovsky on the reconstructed one, with Tafel steps largely potential-insensitive. Together, these results connect atomistic energetics to steady-state polarisation behaviour, clarifying how bias and hydrogen configurations shape HER activity. The trends further suggest that modifying the local electronic environment through electropositive species may reduce critical barriers and enhance performance.

---

1 Introduction

Hydrogen is recognised as a sustainable energy carrier and a key foundation of the chemical value chain, underpinning numerous industrial products. The rise in energy demand and greenhouse gas emissions underscores the urgency of clean alternatives.^1–5 Green hydrogen, produced through renewable-powered water electrolysis, is a promising solution,^6–9 yet the process remains hampered by high overpotentials (η). Overcoming this limitation requires advances in cell efficiency and the discovery of abundant and stable electrocatalysts, both of which are vital to achieve a sustainable hydrogen transition.^7,9,10 Precious metals such as Pt are efficient hydrogen evolution reaction (HER) catalysts, but their scarcity and cost hinder large-scale deployment. Transition metal phosphides (TMPs), particularly Ni₂P, have therefore emerged as promising alternatives, offering abundance, stability, and well-documented catalytic activity.^6,10–30

The catalytic activity of Ni₂P has been demonstrated through predictions of density functional theory (DFT) and experimental validation, with subsequent research directed toward catalyst engineering.^{17,18,25,30–36} Mechanistic studies revealed hydrogen adsorption at specific nickel and phosphorus sites, where phosphorus plays a critical role,^{13,26,37–40} yet systematic computational analyses of HER transition states and activation energies on Ni₂P remain scarce.

The computational hydrogen electrode (CHE) has become a widely used framework, successfully reproducing experimental trends for HER and oxygen evolution reaction (OER).^41,42 However, by treating the potential as a uniform energy shift and neglecting explicit solvent effects, CHE cannot fully capture the electrochemical complexity. Early kinetic studies^43–45 provided valuable information on HER mechanisms but did not incorporate potential-dependent barrier changes or site heterogeneity, reinforcing the need for models that connect atomistic energetics with experimentally measurable responses.

Microkinetic modelling (MKM) has emerged as one such framework, based on pioneering work in heterogeneous catalysis.^46–48 Other complementary approaches, such as kinetic Monte Carlo (kMC), have also been developed for reaction kinetics, providing a stochastic route to capture dynamic processes.⁴⁹ Unlike Butler–Volmer or Marcus theory, MKM, in combination with grand-canonical transition state theory (GC-TST),^50–52 does not impose symmetry parameters and allows rate-determining steps to emerge naturally.^53,54 This makes MKM particularly suited for linking atomistic energetics to potential-dependent activity predictions.

Despite extensive work on transition metals, HER studies on Ni₂P have mainly emphasised the Volmer and Tafel steps.^17,40 Evidence comparing Heyrovsky and Tafel steps is limited and often relies on gas-phase adsorption studies that neglect solvent effects, reducing kinetic precision. Understanding how solvent and potential explicitly affect hydrogen adsorption and transition-state energetics on Ni₂P terminations, and systematically connecting these insights to electrochemical behaviour through MKM, is still lacking.^55–61

In this work, we address these challenges by investigating the HER on two surface terminations of Ni₂P(001) using DFT with constant-potential transition states and a hybrid solvent description. Extensive work on HER-MoS₂^52,60 employs similar methodologies for mechanistic insights. Secondly, microkinetic modelling is applied to link elementary steps to polarisation behaviour, providing a mechanistic framework that bridges atomic-scale theory with electrochemical responses and informs rational catalyst design.^60–62 The present study focuses exclusively on potentials within the stability window of Ni₂P, as the catalyst is reported to be unstable below −0.6 V and above −0.1 V at 0 pH.^13,63–65 As one of the main findings, our results show that the phosphorus reconstruction of Ni₂P(001) leads to a substantial reduction of the HER overpotential, yielding values in excellent agreement with experiments.

The remainder of the paper is organised as follows. Section 2 outlines the theoretical and computational framework, covering the hydrogen evolution pathway, solvent model, constant potential approach, microkinetic formulation, and computational details. Section 3 presents the results, from gas-phase reference calculations to reaction energetics on pristine and reconstructed surfaces, and concludes with the catalytic response obtained from microkinetic modelling. Section 4 summarises the key findings.

2 Methods

2.1 Computational details

The HER under electrochemical conditions was investigated within a framework that combines density functional theory (DFT) energetics, implicit solvation, and constant-potential grand canonical modelling. This approach captures the electrochemical response by linking potential-dependent transition state theory to steady-state microkinetic simulations, providing a consistent treatment from the electronic structure level through to macroscopic kinetics.

All calculations were performed using the Vienna ab initio Simulation Package^66–68 with the PBE functional within the generalised gradient approximation⁶⁹ and projector augmented-wave potentials.⁷⁰ A plane-wave basis with a kinetic energy cutoff of 500 eV was used throughout. Partial occupancies near the Fermi level were treated with a Gaussian smearing of 0.05–0.1 eV depending on the environment. A 6 × 6 × 1 Γ-centered k-points mesh was used to sample the Brillouin zone for all the calculations. Electronic self-consistency was reached to within 10⁻⁶–10⁻⁷ eV using a combination of the blocked-Davidson and RMM-DIIS algorithms.^71–73 Structural optimisations were performed with the fast inertial relaxation engine (FIRE) algorithm,^74,75 achieving force convergence below 0.025 eV Å⁻¹ in all cases. For the gas-phase density of states (DOS) calculations, the only changes from the above parameters were the inclusion of dipole corrections along the z-direction for the total energy and forces, and a denser 11 × 11 × 1 Γ-centred k-point mesh.

The role of dispersion interactions was evaluated, since hydrogen adsorption free energies are central to the mechanistic analysis. Earlier work on Ni₂P showed that adding a D3 correction to PBE leads to an almost constant downward shift of approximately 0.1 eV in the differential Gibbs free energy across hydrogen coverage.^76,77 This shift reflects an overbinding tendency in the PBE+D3 description and aligns with observations for metallic systems, where PBE hydrogen adsorption energies agree more closely with experiment than PBE-D3 or RPBE.^78,79 Turning to water–metal interactions, dispersion corrections are generally regarded as beneficial. However, PBE+D3 tends to overestimate adsorption energies relative to PBE.⁸⁰ RPBE+D3 improves the accuracy of water adsorption and gives hydrogen adsorption energetics that are essentially equivalent to PBE.⁷⁸ These trends show that the choice of dispersion correction and functional involves a compromise between water adsorption and hydrogen adsorption accuracy. Because this work depends primarily on hydrogen adsorption energetics, PBE without D3 was adopted for all reaction calculations. A detailed comparison of these deviations lies beyond the scope of this study.

The catalyst slabs were generated as 2 × 2 surface supercells derived from the converged bulk structure reported in our previous work,⁷⁶ which is consistent with earlier studies.^81–83 For electrochemical calculations, the slabs consisted of five atomic layers separated by a 14 Å solvent region, including the explicit water cluster. The three bottom layers were fixed at their bulk positions, while the two top layers and adsorbates were fully optimised. The DOS calculations were performed using seven-layer slabs with the three bottom layers fixed and a 20 Å vacuum separation. When frozen-geometry calculations were required, all atomic positions were kept fixed. Spin-polarised calculations confirmed that the reference slabs and adsorption structures were non-magnetic at the considered coverages.

The transition states were located using CI-NEB calculations^84–86 with a variable number of intermediate images, followed by refinement with the dimer method.^87–89 Both approaches were employed as implemented in the VTST code suite.⁹⁰ Neutral, optimised transition state geometries were subsequently frozen, and static calculations with fractional electron counts were used to obtain grand canonical electronic energies. Electrolyte effects were included using the VASPsol implicit solvent model^91–93 based on the joint-DFT (jDFT)⁹⁴ with a dielectric constant of 78.4 and a Debye screening length of 3.0 Å. The resulting potential-dependent free energies of initial, final, and transition states were used in the microkinetic modelling of HER.

2.2 Hydrogen evolution

In acidic media, the hydrogen evolution reaction (HER) proceeds through the Volmer, Heyrovsky, and Tafel steps:

H⁺ + * + e⁻ → H* (Volmer) (1)

H* + H⁺ + e⁻ → H₂ + * (Heyrovsky) (2)

2H* → H₂ + 2* (Tafel) (3)

where * denotes a vacant surface site. The overall process was analysed through two competing pathways: (i) sequential Volmer steps followed by a Tafel step (denoted VVT), and (ii) a Volmer step coupled with a Heyrovsky step (denoted VH).

Hydrogen adsorption, reflecting surface affinity, is central to defining reaction pathways. Coverage-dependent differential Gibbs free energies were established in our previous work,⁷⁶ which identified the saturation hydrogen coverages for both pristine and phosphorus-reconstructed Ni₃P₂-terminated Ni₂P(001) surfaces. These results form the thermodynamic basis for the present investigation.

2.3 Solvent model

The solvent environment was described using a hybrid approach that combines implicit solvation with explicit water clusters, thereby accounting for both long-range dielectric screening and short-range hydrogen-bonding effects at the electrode–electrolyte interface. Explicit clusters of either H₉O₄⁺ (Eigen cation) or H₈O₄ were introduced above the catalyst surface. The Eigen cation was employed for the Volmer and Heyrovsky steps, which involve proton transfer from the solvent, while the neutral H₈O₄ cluster was used for the Tafel step.

The implicit contribution was treated within a polarizable continuum model,^92,93,95 where the relative permittivity is expressed as a functional of the solute electron density:

ε(n(r⃑)) = 1 + (ε_b − 1)S(n(r⃑)), (4)

with S(n) defining the smooth transition between solute, cavity, and solvent regions (see SI). This formulation introduces an additional solvent potential V_solv into the Kohn–Sham Hamiltonian and leads to a linearised Poisson–Boltzmann (LPB) description of the electrostatic potential in solution. Further details of the implicit solvent description are provided in the SI.

2.4 Constant potential

Electrochemical reactions occur at constant applied potentials, whereas periodic DFT calculations are performed at constant charge. To reconcile this, a constant-potential framework is adopted^52,59,60,62 in terms of the grand canonical energy,

Ω(U) = E^DFT_elec(U) − δq_elec(U)U − δq_elec(U)V_solv, (5)

where E^DFT_elec is the DFT electronic energy, δq_elec(U) is the electrode charge at potential U and V_solv is the potential in the bulk solvent. The final term in eqn (5) is specific to the VASPSol implementation.^91–93 U and V_solv are relative to the vacuum scale. The applied potential U is defined as U = (−ε_f/|e|) − V_solv. It is computed by varying the electron number δn_e in the simulation cell and adding the bulk solvent potential to the Fermi energy. The potential U is referenced to that of the standard hydrogen electrode (SHE) in vacuum scale. The remainder of the present work uses potential U vs. SHE. This formalism enables the direct mapping of reaction energetics (see Fig. 1) as a function of the applied potential. Details are provided in the SI (Section S1.2).

Fig. 1 Grand canonical energies (Ω) of initial, final, and transition states as a function of applied potential. Activation (reaction) energies at the potential of zero charge (PZC) denoted by Ω^PZC_a(r).

To limit the computational cost of geometry optimisations under this scheme, the transition state (TS) structures were optimised once at the potential of zero charge (PZC) and subsequently kept fixed. For each applied potential, the electron number was varied on this frozen geometry until electronic self-consistency was reached. Thus, the TS structure at all potentials corresponds to that optimised at the PZC. In contrast, all initial and final states were fully geometry optimised at each potential. The validity of the frozen-TS approximation was assessed through benchmark calculations on selected TS searches, with the results presented in the SI (Fig. S1). Furthermore, the impact of dipole interactions associated with asymmetric slab models was benchmarked against symmetric slab models. The Volmer step, taken as a representative charge-transfer reaction with significant dipole changes, shows activation-energy variations of at most 0.06 eV (see SI Sections S1.4.2 and S1.4.3).

The backward Heyrovsky and Tafel steps were neglected to reflect low H₂ conditions, while a 0.2 eV minimum was enforced for the Volmer and Heyrovsky barriers to account for proton transfer and solvent reorganisation at a mean-field level.^62,96 Proton-coupled electron transfer (PCET) at electrified interfaces involves solvent reorganisation and proton motion, which are not explicitly resolved in the present modelling framework. Accordingly, explicit solvent dynamics and proton transport kinetics, as well as more advanced PCET descriptions, lie beyond the scope of the present study and is investigated in detail elsewhere.^97–102

2.5 Microkinetic model

In the microkinetic framework, each hydrogen coverage is treated as a distinct state of the system, with electrochemical reactions enabling transitions between states. The kinetics are described by the GC-TST,⁵¹ where the transition rate constant from states j to i is given byHere, Ω^ij_a denotes the grand-canonical activation barrier. The transmission coefficient was set to unity, consistent with electronically adiabatic reactions at constant potential.^51,62

The steady-state solution of the master equation,

yields the probability of each hydrogen coverage, from which elementary step rates and the net current density are obtained. The current density is expressed aswhere r_v/h⁺ and r_v/h⁻ are the forward and backward Volmer/Heyrovsky rates, e is the elementary charge, and A the surface area per adsorption site. Four and six adsorption sites were used for the pristine and reconstructed surfaces, respectively, corresponding to a simulation cell area of 119.29 × 10⁻¹⁶ cm². This framework enables the construction of theoretical polarisation curves and Tafel plots for both terminations. See SI for details.

3 Results and discussion

The presentation of results builds on the theoretical framework and computational methods, beginning with gas phase modelling (Section 3.1), advancing to reaction energetics on pristine and reconstructed surfaces (Sections 3.2.1 and 3.2.2), and concluding with an assessment of the catalytic response through microkinetic modelling (Section 3.3).

3.1 Hydrogen adsorption

The structural optimisation resulted in a bulk structure of Ni₂P with a hexagonal lattice and lattice parameters a = b = 5.871 Å and c = 3.381 Å.⁷⁶ Following the trend in formation energy as reported,⁷⁶ surface terminations Ni₃P₂ and Ni₃P₂ + 4P for Ni₂P(001) were considered for adsorption and reaction calculations (see Fig. 2). This choice of surface termination is also supported by experiments.^103,104 Hydrogen adsorption was investigated using a 2 × 2 surface with a 7-layered slab, while reaction studies used a 5-layered slab. In both cases, the bottom three layers were fixed.

Fig. 2 Side views of (a) Ni₃P₂ and (d) Ni₃P₂ + 4P surface terminations of Ni₂P(001). Red, blue and white spheres denote Ni, P and H atoms, respectively. Hydrogen-saturated surfaces are shown for (c) Ni₃P₂ (4H) and (d) Ni₃P₂ + 4P (6H) in a 2 × 2 unit cell. Hydrogen (4 atoms) occupies the available hollow sites on the pristine surface. For the reconstructed surface, the lowest energy configuration for 6H atoms corresponds to a structure where 3 out of 4 external P atoms bind 2H, each.⁷⁶

Previously, our systematic analysis of the differential Gibbs free energy⁷⁶ suggested a coverage of 4 and 6 hydrogen atoms as saturation adsorption on pristine and reconstructed 2 × 2 surfaces, respectively (see Fig. 2). The pristine surface is saturated with highly exothermic hydrogen adsorption, whereas the reconstructed surface is saturated close to thermoneutral hydrogen adsorption.

In addition to thermodynamic trends of adsorption, the density of states (DOS) was examined to assess how surface reconstructions modify the electronic structure. Bulk Ni₂P and the total DOS of symmetric 15- and 7-layer Ni₃P₂-terminated slabs confirm metallic character resembling the bulk (Fig. S5a and b), providing a benchmark for surface analyses.

For the 7-layer Ni₃P₂ termination, the pristine surface largely retains metallicity, with Ni 3d orbitals having the largest contribution near the Fermi level (see Fig. 3a, S5d). As the surface reconstructs, the predominant Ni 3d peak around −1.5 eV is reduced and shifted toward deeper valence bands, and the Fermi-level DOS is nearly halved, though metallic character is still preserved (see Fig. 3b).¹⁰⁵

Fig. 3 Atom-resolved density of states (DOS) for (a) Ni₃P₂, and (b) Ni₃P₂ + 4P terminations of Ni₂P(001) surface.

Hydrogen adsorption sharpens this contrast. On pristine surfaces, it further depletes the Fermi-level DOS, while on reconstructed surfaces the profile remains essentially unchanged (Fig. S4). Metallic conductivity therefore persists across all relevant coverages, ensuring efficient charge transfer. Combined with the availability of P adatom sites, this offers particularly favourable conditions for HER.

3.2 Reaction calculations

The HER pathways on the selected surface terminations are analysed through the VVT and VH mechanisms, denoted PX_vvt/vh and RX_vvt/vh for the pristine and reconstructed surface, respectively, where X denotes a number. For a detailed discussion of the reaction mechanisms, the initial, transition, and final states of the Volmer, Tafel, and Heyrovsky steps are denoted as IS_V/T/HX, TS_V/T/HX, and FS_V/T/HX, respectively. These reaction cycles, intermediate steps, and energetics are summarised in Tables 1 and 2. Given the consistent trends in activation and reaction energies in negative applied potentials, the discussion focuses on a representative potential of −0.4 V. This value provides a clear basis for analysing mechanistic responses under cathodic bias, where most reaction cycles are exothermic. The initial state IS_V1 differs in energy between the two voltages shown in the figures, producing a constant shift in the energy profiles. For clarity, all energy curves are referenced to 0 eV at the starting point.

Table 1 Activation (Ω^cycle_a) and reaction (Ω^cycle_r) free energies at −0.4 V (0 V) for complete HER cycles via Volmer–Volmer–Tafel (VVT) and Volmer–Heyrovsky (VH) pathways on pristine (PX_vvt/vh) and reconstructed (RX_vvt/vh) surfaces. The values of the best pathways are denoted in boldface for each surface

Cycle	Surface coverage transition	Ω^cyclea (eV)	Ω^cycler (eV)
P1vvt	3H(0) → 4H(0) → 5H(0) → 3H(0)	0.37 (0.53)	−0.21 (0.49)
P2vvt	4H(0) → 5H(I) → 6H(II) → 4H(I)	0.57 (1.19)	0.48 (1.14)
P3vvt	4H(0) → 5H(I) → 6H(I) → 4H(0)	0.52 (0.96)	−0.24 (0.49)
P4vh	4H(0) → 5H(I) → 4H(0)	0.52 (0.77)	−0.24 (0.47)
P5vh	4H(0) → 5H(0) → 4H(0)	0.52 (0.78)	−0.2 (0.54)
R6vvt	6H(0) → 7H(0) → 8H(0) → 6H(0)	1.02 (1.63)	−0.15 (0.54)
R7vh	6H(0) → 7H(I) → 6H(0)	0.65 (0.88)	−0.14 (0.47)
R8vh	6H(0) → 7H(0) → 6H(0)	0.23 (0.56)	−0.13 (0.55)

---

Table 2 Activation (Ω_a) and reaction (Ω_r) energies at −0.4 V (0 V) for Volmer, Tafel, and Heyrovsky steps at various coverages on pristine (Ni₃P₂) and reconstructed (Ni₃P₂ + 4P) surfaces. The energetically most favourable values are denoted in boldface for each reaction step

		Configurations		GCE (eV)
		Initial	Final	Ωa	Ωr
Pristine	V	3H(0)	4H(0)	0.37 (0.53)	−0.42 (−0.09)
		4H(0)	5H(I)	0.52 (0.72)	0.25 (0.57)
		4H(0)	5H(0)	0.52 (0.72)	0.27 (0.58)
		5H(I)	6H(I)	0.67 (0.88)	0.41 (0.73)
		5H(I)	6H(II)	0.52 (0.75)	0.24 (0.59)
	T	5H(0)	3H(0)	0.59 (0.52)	0.51 (0.48)
		6H(I)	4H(0)	0.24 (0.15)	−0.39 (−0.32)
		6H(II)	4H(I)	0.55 (0.49)	0.46 (0.44)
	H	5H(I)	4H(0)	0.26 (0.65)	0.03 (0.43)
		5H(0)	4H(0)	0.24 (0.61)	0.04 (0.44)
Reconstructed	V	6H(0)	7H(I)	0.58 (0.72)	0.45 (0.61)
		6H(0)	7H(0)	0.19 (0.28)	0.03 (0.23)
		7H(0)	8H(0)	0.54 (0.76)	0.13 (0.47)
	T	8H(0)	6H(0)	1.39 (1.41)	0.22 (0.32)
	H	7H(I)	6H(0)	0.35 (0.69)	−0.10 (0.28)
		7H(0)	6H(0)	0.34 (0.77)	0.32 (0.77)

---

3.2.1 Pristine surface – Ni₃P₂. Among the VVT pathways on the pristine surface, P1_vvt (see Fig. 4) is energetically the most favourable (see Table 1). Within this cycle, one of the hollow sites is initially vacant (3H configuration), being filled by the first Volmer step, followed by another one on a nearby bridge site. After the Tafel step (5H), H₂ desorption follows and the system returns to its original 3H coverage (see Fig. 4e). Although energetically the most favourable, this pathway involves events that proceed thermodynamically uphill, even at an applied potential of −0.4 V.

Fig. 4 Hydrogen evolution via the P1_vvt pathway on the Ni₃P₂ surface. The orange spheres denote oxygen atoms, and the remaining colours are defined in Fig. 2. IS, TS, and FS are shown for the Volmer steps from (a) 3 to 4 and (b) 4 to 5H coverages, and for (c) H₂ formation and (d) desorption via the Tafel mechanism, and (e) a consolidated illustration. The explicit water cluster is surrounded by implicit solvent.

As shown in Fig. 5, the free-energy landscape decreases until the proton-transfer step and then ascends, indicating that H₂ desorption involves a higher barrier than its formation. Under negative bias, IS_V2 is markedly stabilised relative to IS_V1, producing a deep potential well, an effect attributed to the strong electrostatic interaction between the saturated 4H surface and the Eigen cation.

Fig. 5 Grand-canonical energy landscapes of P1_vvt for two applied potentials 0 V (black) and −0.4 V (red). The evolutionary stages of H are represented by H⁺ (dissolved proton), H* (adsorbed H), H₂* (adsorbed H₂), and H₂ (desorbed H₂). TS_v/t/H⁺ denotes the transition states of the Volmer, Tafel, and proton-transfer steps from the solvent to the interfacial water cluster. The vertical arrows indicate activation/reaction energy Ω_a/r of the corresponding cycle.

The P3_vvt pathway that begins with a saturated 4H coverage (see Fig. S11d) does not exhibit a uniformly downward free-energy landscape, although stabilisation of the 5H intermediate via Eigen cation interaction improves the overall thermodynamics. Along the reaction coordinate, the Tafel step shows a net downward trend with moderate barriers, rendering the 4H surface with desorbed H₂ more stable than the 6H configuration. However, the 5H–6H Volmer step, has a high barrier and remains the rate-limiting step. P2_vvt presents a less demanding 5H–6H Volmer barrier but is overall endothermic (see Fig. S9d), with an uphill Tafel step arising from a difference in the final 4H configuration (bridge site) relative to P3_vvt.

Comparison of P2_vvt and P3_vvt further highlights the influence of surface hydrogen arrangements (see surface coverage transitions in Table 1). In P2_vvt, the second transition state (TS_V2) and the corresponding final state (FS_V2) lie about 0.16 eV lower in energy than in P3_vvt due to alternative 6H arrangements. Hydrogen diffusion between bridge sites can convert a 6H(II) to a 6H(I) configuration, but the associated 0.4 eV energy cost does not offer a significant improvement in the free-energy profile. Although the Tafel steps T_P2 and T_P3 start and end at identical coverages, their thermodynamic character differs: P2_vvt is endothermic, whereas P3_vvt is exothermic (see Tables 1 and 2, Fig. S9d, S11d). These findings emphasise the decisive role of surface hydrogen configurations in the shaping of reaction energetics and the need to identify global minimum coverage states.⁷⁶ Although a systematic search was performed, the vast configurational space remains challenging to explore exhaustively, and future studies can benefit from active-learning approaches to accelerate this process.

Within the VH pathways (see Fig. S12 and S14), P4_vh and P5_vh exhibit nearly identical free-energy profiles, reflecting the comparable adsorption characteristics of the 5H(0) and 5H(I) configurations (see Fig. 6a and b). The Heyrovsky steps are thermodynamically more favourable than the Volmer steps (see Table 2), which remain limiting under moderately negative potentials. These trends suggest that dopants capable of narrowing the Gibbs free-energy gap between saturated and adjacent hydrogen coverages could enhance HER activity. Electropositive dopants, in particular, may reduce barriers for both Volmer and Heyrovsky steps and thereby enhance overall HER performance. They may also benefit the Tafel process, as the mild deterioration of its barriers under negative bias, while largely insensitive, is consistent with the absence of charge transfer (see eqn (3)).

Fig. 6 A consolidated illustration of (a) P4_vh, (b) P5_vh, (c) R7_vh, and (d) R8_vh HER mechanisms.

Taken together, the results for the pristine surface reveal that the HER activity on Ni₃P₂ is primarily constrained by high-barrier Volmer steps and, in certain cycles, by unfavourable Tafel energetics linked to specific hydrogen configurations (see Table 2). Negative potentials preferentially stabilise intermediates and transition states in Heyrovsky steps, indicating that bias sensitivity differs among mechanisms. The Eigen cation exhibits near-degeneracy across multiple spatial configurations within about 0.01 eV, implying that its position does not significantly alter reaction energetics and supporting the use of a fixed water-cluster model. From a design perspective, improving HER performance on Ni₃P₂ will require strategies that lower the transition states of the high-barrier Volmer steps and optimise surface hydrogen configurations.

3.2.2 Reconstructed surface – Ni₃P₂ + 4P. Despite the overall exothermicity of the R6_vvt cycle, which starts from saturated 6H coverage, its activation barriers are prohibitively high and exhibit only weak sensitivity to applied potential (see Fig. S17d). Structural analysis reveals close similarity between the transition states for H₂ recombination and desorption in the Tafel steps (see Fig. S16). Both involve surface diffusion of a PH₂ group to yield an adjacent but sub-optimal 8H configuration. The transition states for forming this configuration and for H₂ release are similar in structure, and both remain energetically costly. Given the combination of high, potential-insensitive barriers and competing low-energy but sub-optimal configurations, R6_vvt is unlikely to contribute to HER, even under moderate negative bias (see Table 1).

R7_vh and R8_vh differ only by the adsorption site of a single hydrogen atom (see Fig. 6c, d, 7, S18). In 7H(0), the additional hydrogen binds to the phosphorus adatom, whereas in 7H(I), it occupies the Ni–P bridge on the pristine layer. The Volmer step that forms 7H(0) is markedly more favourable than that leading to 7H(I) (see Table 2). Consequently, the Volmer step determines the rate in R7_vh, while the Heyrovsky step determines the rate in R8_vh which emerges as the most viable pathway for H₂ formation on the reconstructed surface.

Fig. 7 Hydrogen evolution via the R8_vh pathway on Ni₃P₂ + 4P, showing IS, TS, and FS for (a) the Volmer step from 6 to 7H coverage and (b) the Heyrovsky step from 7 to 6H with H₂ desorption.

Potential-dependent energetics clarify this preference (see Fig. 8, S19d). At 0 V, the Heyrovsky reaction enthalpy (H_R7, H_R8) is endothermic in both configurations, making H₂ release thermodynamically unfavourable. Applying a negative bias renders H_R7 exothermic, enabling the reaction to proceed (see Table 2). This shift indicates that at 0 V, the system favours the 7H(I) state with a solvated proton over the desorbed product (6H + H₂), while at −0.4 V, the negative bias promotes electron transfer from the electrode, thereby stabilising the desorbed H₂ state.

Fig. 8 Grand-canonical energy landscapes of R8_vh for two applied potentials 0 V (black) and −0.4 V (red). The evolutionary stages of H is represented by H⁺ (dissolved proton), H* (adsorbed H), H₂* (adsorbed H₂) and H₂ (desorbed H₂). TS_v/t/H⁺ denotes the transition states of the Volmer, Tafel, and proton-transfer steps from the solvent to the interfacial water cluster.

The reconstructed surface therefore follows the stability order: 7H(0) (Eigen cation) > 8H(0) (neutral water cluster) > 6H(0) + H₂ (neutral water cluster) > 7H(I) (Eigen cation).

Constant-potential energy landscapes indicate that VVT dominates on pristine Ni₃P₂ terminations, whereas VH mechanisms are more favourable on reconstructed surfaces. Under negative potentials, Heyrovsky steps are the most accessible electron-transfer events, followed by Volmer steps, while Tafel steps remain largely insensitive to bias. Although both Volmer and Heyrovsky steps benefit from negative bias, the Volmer step typically determines the overall rate.

Electropositive dopants may enhance HER activity on reconstructed Ni₃P₂ by stabilising transition states for the high-barrier Volmer steps, tuning adsorption site energetics, and improving the local environment for Tafel recombination. However, the investigation of dopant effects is non-trivial and demands careful and detailed analysis. Further, the negligible response of reaction energetics (at fixed coverage) to variations in water cluster geometry supports the use of a simplified solvation model in computational screening, allowing for a broader exploration of dopant and surface-structure effects.

3.3 Microkinetic modelling

The GC-TST provides activation barriers to estimate the energetic feasibility of individual reaction steps, but it does not capture how these steps operate collectively under electrochemical conditions. To address this limitation, we constructed a steady-state microkinetic model (MKM) using GC-TST-derived rate constants (see Section S1.3).

The MKM accounts for the interplay between parallel pathways, coverage-dependent kinetics, and potential effects, enabling quantitative predictions of HER polarisation behaviour and mechanistic trends beyond those obtainable from static energy profiles.

Experimental linear sweep voltammetry studies on Ni₂P report a wide range of η across different surfaces, morphologies, and measurement conditions (Table 3). As most of the values in literature refer to the reversible hydrogen electrode (RHE), we have converted them to the SHE scale for direct comparison with our results, as detailed in Section S1.5.

Table 3 Overpotentials η_j (mV vs. SHE) at various current densities j (mA cm⁻²). Ni₂P (PC) denotes polycrystalline Ni₂P

Surface	Overpotential (mV vs. SHE)
	η10	η20	η100
Ni3P2 (this work)	330	345	370
Ni3P2 + 4P (this work)	120	139	180
Ni2P(001)^18	—	147	214
Ni2P(001)^30	—	—	343 ± 22
Ni2P (111)^12	96	—	∼167
Ni2P (PC)^106	130	—	468

---

In this context, our MKM (see Fig. 9) shows that the phosphorus-reconstructed termination of Ni₂P(001) yields η of 120, 139, and 180 mV at 10, 20, and 100 mA cm⁻², respectively, which are substantially lower than those of the pristine termination (330, 345, and 370 mV). The HER η on the reconstructed surface agrees well with experiment,¹⁸ while that reported in another study³⁰ closely matches the pristine surface. These comparisons suggest that the HER-active surface in ref. 18 is likely reconstructed, whereas that in ref. 30 may correspond to the pristine termination. Although detailed surface-termination data are not available from these works, the observed correlations suggest the active surfaces and their distinct mechanistic behaviour. Despite the Pourbaix diagram²⁶ predicting the reconstructed surface to dominate only at more negative potentials (around −210 mV), our results, together with the experiment, indicate that even a partially reconstructed surface can sustain markedly superior HER kinetics well above this potential.

Fig. 9 Polarisation curves of HER on pristine (Ni₃P₂) and reconstructed (Ni₃P₂ + 4P) surface terminations. P/RX_vh stands for pristine/reconstructed surface Volmer–Heyrovsky cycles. For each surface, the solid line represents the total current density including all explored cycles, while dashed lines correspond to individual cycles. The orange (P4_vh) and purple (P5_vh) curves lie close to the red (Ni₃P₂) curve.

In comparison with literature, the active surface for HER on Ni₂P(001) reported^26,77 corresponds to a hydrogen-adsorbed pristine surface. In contrast, in our case the reaction proceeds via the VH mechanism on the reconstructed surface, consistent with the interpretation of ref. 13. The η obtained in the present work lie within the stability window of Ni₂P,^{13,26,63–65} where phosphorus reconstructions are known to be stable under negative potentials in acidic media.^26,63 Furthermore, the η of polycrystalline and Ni₂P(111) surfaces are comparable to that of the reconstructed Ni₂P(001) termination (see Table 3), suggesting that phosphorus reconstruction could also occur on these facets and contribute to HER activity. Together, these findings emphasise the broader relevance of reconstructed surfaces in governing the catalytic performance of Ni₂P.^25,82

The corresponding Tafel slopes are 42 and 58 mV dec⁻¹ for the pristine and reconstructed surfaces (Fig. S24), respectively, in the higher η range of 167–217 mV. At lower η (42–142 mV), the pristine surface does not show a clear linear relation between η and log|j|, making the extracted slope unreliable, whereas the reconstructed termination follows a well-defined linear trend with a slope of 53 mV dec⁻¹. Popczun et al.¹⁸ reported experimental slopes of about 46 mV dec⁻¹ in the lower range and 81 mV dec⁻¹ in the higher range, closely matching the reconstructed surface at low η but deviating at higher η.

As the η increases, the reconstructed surface shows a modest rise in slope, while the pristine surface decreases sharply, indicating faster HER kinetics on the latter at higher potentials (Fig. S24). The wide-range Tafel behaviour reported by Kong et al.³⁰ is qualitatively reproduced by the pristine surface in the present work, whereas the reconstructed termination exhibits a distinctly different trend. The slope is approximately 30 mV dec⁻¹ lower than the experimental values³⁰ in the 217–292 mV range, suggesting that the MKM model modestly overestimates the pristine-surface HER rate in this region. Overall, the MKM results indicate that surface reconstruction enhances HER kinetics at low η, while the pristine surface becomes more active at higher potentials, underscoring the beneficial role of reconstruction.

On the pristine surface, the most kinetically effective isolated cycles are P4_vh and P5_vh, both of which are predominant across the potential range by the 5H–4H Heyrovsky step. Up to approximately −0.4 V, the 4H–5H forward Volmer and its reverse remain in near equilibrium. Beyond this point, the forward Volmer and Heyrovsky steps accelerate while the reverse Volmer slows, producing a marked rise in current density. This sustained Heyrovsky activity, supported by favourable Volmer kinetics, explains the strong performance of P4_vh and P5_vh in both isolated and global contexts.

The GC-TST barriers suggest that P1_vvt > P3_vvt > P2_vvt in kinetic favourability, yet microkinetic modelling reverses this order to P3_vvt > P2_vvt > P1_vvt. Similar discrepancies between the mechanistic feasibility inferred from activation energies and those derived from kinetic analysis have been reported.⁴⁹ In P1_vvt, a low-barrier 3H–4H Volmer is offset by a sluggish Tafel step, keeping forward and reverse Volmer transitions in almost perfect balance at low and moderate η. Only at high negative voltages does the faster Tafel step disrupt this balance, enabling net forward electron transfer via the Volmer step. Despite the higher Volmer barrier, P3_vvt benefits from a low Tafel barrier that brings the 6H–4H Tafel rate in line with the 5H–6H Volmer rate, allowing net forward transfer at lower η than P2_vvt, where the Tafel step remains rate-limiting.

Global steady-state probabilities and rates differ subtly from isolated behaviour. In the whole network, the occupancies retain the key trends of P4_vh, P5_vh, and P3_vvt, with mutual influence blurring their signatures. As such, probabilities alone cannot unambiguously separate their contributions. However, rate analysis shows that for VVT cycles, network coupling can amplify the advantage of a fast Tafel step (P3_vvt) or further suppress slow-step cycles (P1_vvt). While GC-TST establishes intrinsic barriers, incorporating steady-state microkinetic coupling reveals the kinetic interplay that underpins the experimentally observed current density, with P4_vh and P5_vh emerging as the primary contributors (see Fig. 9).

On the reconstructed surface, the full reaction network analysis shows that the R8_vh cycle is the primary contributor to the total current density, with other cycles becoming appreciable only at much higher η. Yet they contribute to the total polarization curve as the complete network achieves slightly lower η than the best-performing isolated cycle, R8_vh, indicating a synergistic benefit of coupling multiple pathways. This becomes evident for the isolated R7_vh cycle, where the 7th hydrogen adsorbs at a Ni–P bridge site of the pristine-like layer, making it energetically less favourable than R8_vh. However, its Heyrovsky rate constant exceeds that of R8_vh down to about −0.3 V (see Fig. 10). On the other hand, a slow 6H–7H forward Volmer step combined with a fast 7H–6H reverse Volmer step maintains a high steady-state population of 6H and a low population of 7H. As a result, the 7H–6H Heyrovsky step contributes little to the current, especially at lower voltages. In the complete network, increased 7H probabilities activate this step down to about −0.3 V, enhancing the rate of H₂ conversion and contributing to lower η compared to isolated R8_vh.

Fig. 10 Comparison of the R7_vh and R8_vh rate constants. The inset shows a magnified view highlighting the crossover between the two cycles.

The global steady-state distribution closely resembles that of R6_vvt (see Fig. S17a and S23c), with notable differences for the 6H and 7H states. In the branch R8_vh –R6_vvt of the combined network, the highly accessible R8_vh pathway slightly stabilises 6H at higher voltages compared to isolated R6_vvt. At lower negative voltages, 6H is the most stable and 8H the least stable, but this ordering reverses as the potential becomes more negative. Under these conditions, the Heyrovsky rate in the complete network slows relative to R8_vh due to reduced 6H and 7H probabilities. A similar pattern occurs in R6_vvt, where competing transitions, 7H–8H forward Volmer and 8H–7H reverse Volmer, become kinetically faster with increasing η. These processes are also active in the combined network, producing a composite rate profile with features of R6_vvt and R8_vh.

This kinetic interplay influences current density primarily at higher voltages. However, the η required to reach 10 mA cm⁻² and 100 mA cm⁻² remain well below this range. The analysis thus highlights how coupling reaction cycles through contributions from R7_vh at moderate voltages and kinetically accessible alternative transitions at higher voltages can slightly outperform the best isolated cycle. A more comprehensive understanding will require mapping a broader set of pathways involving diverse hydrogen coverage configurations.

4 Conclusions

In order to elucidate hydrogen evolution pathways on Ni₂P(001), the present theoretical investigation employs DFT, a hybrid solvation model, grand-canonical transition-state theory and microkinetic modelling, revealing how surface reconstruction and operating conditions shape distinct reaction mechanisms.

The GC-TST analysis shows that the pristine Ni₂P(001) surface favours a Volmer–Volmer–Tafel sequence, limited by high Volmer barriers and sensitive to adsorption coverage configurations, while reconstructed Ni₃P₂ terminations create chemical environments where Volmer–Heyrovsky becomes the most accessible route, with Heyrovsky limiting HER. These contrasts highlight the mechanistic diversity introduced by reconstructions and the importance of mapping the entire energy landscape for each pathway.

Microkinetic modelling reveals that, under steady-state conditions, both pristine and phosphorus-reconstructed Ni₂P(001) surfaces are governed by Volmer–Heyrovsky cycles, in contrast to the GC-TST prediction for pristine terminations. This shows that TST is necessary for dissecting the elementary-step energetics but insufficient on its own to predict polarisation behaviour without steady-state modelling. Additionally, the steady-state hydrogen evolution is limited by the Volmer step on the pristine surface and the Heyrovsky step on the reconstructed surface. The polarisation curves and low η obtained for the reconstructed surface (120–180 mV at j = 10–100 mA cm⁻²) are in good agreement with the available experimental data and highlight the good HER performance of Ni₂P. On the other hand, the (ideal) pristine surface displays substantially higher η, demonstrating the importance of surface morphology on electrochemical performance.

Finally, constant-potential investigations on HER suggest that informed doping of Ni₂P can enhance performance. Electropositive dopants may lower Volmer barriers and optimise the Heyrovsky step, while electronegative dopants may benefit the Tafel step. For realistic polycrystalline catalysts, predictive HER modelling requires integration of multiple facets, surface reconstructions, and full hydrogen coverage maps, ideally with solvent and constant-potential effects. Active-learning-driven coverage exploration has great potential to help overcome key computational bottlenecks, supporting more reliable and broadly applicable simulation-guided catalyst design.

Conflicts of interest

There are no conflicts to declare.

Data availability

Supplementary information is available. See DOI: https://doi.org/10.1039/d5cp04301h.

The constant-potential DFT and microkinetic modelling data used to produce this work are available at Zenodo.¹⁰⁷ The codes used to analyse the data are available at https://github.com/syamsadangit/scripts.

Acknowledgements

The simulations were carried out using the supercomputing resources provided by Sigma2 – the National Infrastructure for High-Performance Computing and Data Storage in Norway (project no. NN9497K) and CSC – IT Centre for Science, Finland (project no. tty1701).

References

1. O. O. Yolcan, Innovation Green Dev., 2023, 2, 100070.
2. Q. Wang, T. Yang and R. Li, Environ. Res., 2023, 216, 114575.
3. M. Mahdavi and D. Vera, Int. J. Hydrogen Energy, 2023, 48(88), 34575–34598.
4. F. Chien, L. Huang and W. Zhao, J. Innovation Knowl., 2023, 8, 100298.
5. C. M. S. Kumar, S. Singh, M. K. Gupta, Y. M. Nimdeo, R. Raushan, A. V. Deorankar, T. A. Kumar, P. K. Rout, C. Chanotiya and V. D. Pakhale, et al., Sustainable Energy Technol. Assess., 2023, 55, 102905.
6. Z. Pu, T. Liu, I. S. Amiinu, R. Cheng, P. Wang, C. Zhang, P. Ji, W. Hu, J. Liu and S. Mu, Adv. Funct. Mater., 2020, 30, 2004009.
7. V. Maestre, A. Ortiz and I. Ortiz, Renewable Sustainable Energy Rev., 2021, 152, 111628.
8. R. Lowe and P. Drummond, Renewable Sustainable Energy Rev., 2022, 153, 111720.
9. V. Panchenko, Y. V. Daus, A. Kovalev, I. Yudaev and Y. V. Litti, Int. J. Hydrogen Energy, 2023, 48, 4551–4571.
10. Q. Zhou, Z. Shen, C. Zhu, J. Li, Z. Ding, P. Wang, F. Pan, Z. Zhang, H. Ma and S. Wang, et al., Adv. Mater., 2018, 30, 1800140.
11. S. Fujii, S. Ishida and S. Asano, J. Phys. F: Met. Phys., 1988, 18, 971.
12. L. Xiong, B. Wang, H. Cai, H. Hao, J. Li, T. Yang and S. Yang, Appl. Catal., B, 2021, 295, 120283.
13. R. B. Wexler, J. M. P. Martirez and A. M. Rappe, ACS Catal., 2017, 7, 7718–7725.
14. P. Ji, H. Jin, H. Xia, X. Luo, J. Zhu, Z. Pu and S. Mu, ACS Appl. Mater. Interfaces, 2019, 12, 727–733.
15. W. Gao, M. Yan, H.-Y. Cheung, Z. Xia, X. Zhou, Y. Qin, C.-Y. Wong, J. C. Ho, C.-R. Chang and Y. Qu, Nano Energy, 2017, 38, 290–296.
16. R. Ge, J. Huo, T. Liao, Y. Liu, M. Zhu, Y. Li, J. Zhang and W. Li, Appl. Catal., B, 2020, 260, 118196.
17. P. Liu and J. A. Rodriguez, J. Am. Chem. Soc., 2005, 127, 14871–14878.
18. E. J. Popczun, J. R. McKone, C. G. Read, A. J. Biacchi, A. M. Wiltrout, N. S. Lewis and R. E. Schaak, J. Am. Chem. Soc., 2013, 135, 9267–9270.
19. Y. Xu, R. Wu, J. Zhang, Y. Shi and B. Zhang, Chem. Commun., 2013, 49, 6656–6658.
20. J. Tian, Q. Liu, A. M. Asiri and X. Sun, J. Am. Chem. Soc., 2014, 136, 7587–7590.
21. Q. Liu, J. Tian, W. Cui, P. Jiang, N. Cheng, A. M. Asiri and X. Sun, Angew. Chem., Int. Ed., 2014, 53, 6710–6714.
22. Z. Pu, Q. Liu, A. M. Asiri and X. Sun, ACS Appl. Mater. Interfaces, 2014, 6, 21874–21879.
23. J. Tian, Q. Liu, N. Cheng, A. M. Asiri and X. Sun, Angew. Chem., 2014, 126, 9731–9735.
24. J. Kibsgaard, C. Tsai, K. Chan, J. D. Benck, J. K. Nørskov, F. Abild-Pedersen and T. F. Jaramillo, Energy Environ. Sci., 2015, 8, 3022–3029.
25. A. B. Laursen, R. B. Wexler, M. J. Whitaker, E. J. Izett, K. U. Calvinho, S. Hwang, R. Rucker, H. Wang, J. Li and E. Garfunkel, et al., ACS Catal., 2018, 8, 4408–4419.
26. S. Banerjee, A. Kakekhani, R. B. Wexler and A. M. Rappe, ACS Catal., 2023, 13, 4611–4621.
27. S. Anantharaj, S. R. Ede, K. Sakthikumar, K. Karthick, S. Mishra and S. Kundu, ACS Catal., 2016, 6, 8069–8097.
28. Y. Pan, Y. Liu, J. Zhao, K. Yang, J. Liang, D. Liu, W. Hu, D. Liu, Y. Liu and C. Liu, J. Mater. Chem. A, 2015, 3, 1656–1665.
29. S. Cao, Y. Chen, C.-J. Wang, P. He and W.-F. Fu, Chem. Commun., 2014, 50, 10427–10429.
30. S. Kong, D. Raturi, B. Owens-Baird, W. Zheng, Y. V. Kolen’ko, P. Singh, D. D. Johnson and K. Kovnir, ACS Catal., 2025, 15(21), 18723–18737.
31. L. Feng, H. Vrubel, M. Bensimon and X. Hu, Phys. Chem. Chem. Phys., 2014, 16, 5917–5921.
32. C.-C. Weng, J.-T. Ren and Z.-Y. Yuan, ChemSusChem, 2020, 13, 3357–3375.
33. M. H. Hansen, L.-A. Stern, L. Feng, J. Rossmeisl and X. Hu, Phys. Chem. Chem. Phys., 2015, 17, 10823–10829.
34. Z. Wang, S. Wang, L. Ma, Y. Guo, J. Sun, N. Zhang and R. Jiang, Small, 2021, 17, 2006770.
35. C. Liu, T. Gong, J. Zhang, X. Zheng, J. Mao, H. Liu, Y. Li and Q. Hao, Appl. Catal., B, 2020, 262, 118245.
36. L. Chen and X.-K. Wei, Energies, 2024, 17, 2294.
37. H. Ariga, M. Kawashima, S. Takakusagi and K. Asakura, Chem. Lett., 2013, 42, 1481–1483.
38. L. Partanen, M. Hakala and K. Laasonen, Phys. Chem. Chem. Phys., 2019, 21, 184–191.
39. Q. Yuan, H. Ariga and K. Asakura, Top. Catal., 2015, 58, 194–200.
40. J.-S. Moon, J.-H. Jang, E.-G. Kim, Y.-H. Chung, S. J. Yoo and Y.-K. Lee, J. Catal., 2015, 326, 92–99.
41. J. K. Nørskov, J. Rossmeisl, A. Logadottir, L. Lindqvist, J. R. Kitchin, T. Bligaard and H. Jonsson, J. Phys. Chem. B, 2004, 108, 17886–17892.
42. J. K. Nørskov, T. Bligaard, A. Logadottir, J. Kitchin, J. G. Chen, S. Pandelov and U. Stimming, J. Electrochem. Soc., 2005, 152, J23.
43. J. Bockris and E. Potter, J. Electrochem. Soc., 1952, 99, 169.
44. R. Parsons, Trans. Faraday Soc., 1958, 54, 1053–1063.
45. R. Parsons, Trans. Faraday Soc., 1960, 56, 1340–1350.
46. M. Boudart, AIChE J., 1956, 2, 62–64.
47. H. Lynggaard, A. Andreasen, C. Stegelmann and P. Stoltze, Prog. Surf. Sci., 2004, 77, 71–137.
48. K. Honkala, A. Hellman, I. Remediakis, A. Logadottir, A. Carlsson, S. Dahl, C. H. Christensen and J. K. Nørskov, Science, 2005, 307, 555–558.
49. K. Reuter and M. Scheffler, Phys. Rev. B:Condens. Matter Mater. Phys., 2006, 73, 045433.
50. M. M. Melander, Curr. Opin. Electrochem., 2021, 29, 100749.
51. M. M. Melander, J. Electrochem. Soc., 2020, 167, 116518.
52. N. Abidi, A. Bonduelle-Skrzypczak and S. N. Steinmann, Int. J. Hydrogen Energy, 2023, 48, 8478–8488.
53. A. Baz, S. T. Dix, A. Holewinski and S. Linic, J. Catal., 2021, 404, 864–872.
54. A. J. Medford, M. R. Kunz, S. M. Ewing, T. Borders and R. Fushimi, ACS Catal., 2018, 8, 7403–7429.
55. A. J. Garza, A. T. Bell and M. Head-Gordon, ACS Catal., 2018, 8, 1490–1499.
56. A. J. Garza, A. T. Bell and M. Head-Gordon, J. Phys. Chem. Lett., 2018, 9, 601–606.
57. M. R. Singh, J. D. Goodpaster, A. Z. Weber, M. Head-Gordon and A. T. Bell, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, E8812–E8821.
58. F. Yu, H. Zhou, Y. Huang, J. Sun, F. Qin, J. Bao, W. A. Goddard III, S. Chen and Z. Ren, Nat. Commun., 2018, 9, 2551.
59. M. Van den Bossche, E. Skúlason, C. Rose-Petruck and H. Jónsson, J. Phys. Chem. C, 2019, 123, 4116–4124.
60. S. Ø. Hanslin, H. Jónsson and J. Akola, Phys. Chem. Chem. Phys., 2023, 25, 15162–15172.
61. J. D. Goodpaster, A. T. Bell and M. Head-Gordon, J. Phys. Chem. Lett., 2016, 7, 1471–1477.
62. S. Ø. Hanslin, H. Jónsson and J. Akola, ChemPhysChem, 2024, 25, e202400349.
63. A. R. Kucernak and V. N. N. Sundaram, J. Mater. Chem. A, 2014, 2, 17435–17445.
64. A. Parra-Puerto, K. L. Ng, K. Fahy, A. E. Goode, M. P. Ryan and A. Kucernak, ACS Catal., 2019, 9, 11515–11529.
65. R. A. Rivera-Maldonado, A. J. Gironda, A. Varughese, D.-Y. Kuo, H. Nguyen, D. Dean-Hill, J. E. Abramson, G. T. Seidler and B. M. Cossairt, J. Phys. Chem. C, 2025, 129, 1165–1172.
66. G. Kresse and J. Hafner, Phys. Rev. B:Condens. Matter Mater. Phys., 1994, 49, 14251–14269.
67. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
68. G. Kresse and J. Furthmüller, Phys. Rev. B:Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
69. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
70. G. Kresse and D. Joubert, Phys. Rev. B:Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
71. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
72. G. Kresse and J. Furthmüller, Phys. Rev. B:Condens. Matter Mater. Phys., 1996, 54, 11169.
73. P. Pulay, Chem. Phys. Lett., 1980, 73, 393–398.
74. E. Bitzek, P. Koskinen, F. Gähler, M. Moseler and P. Gumbsch, Phys. Rev. Lett., 2006, 97, 170201.
75. D. Sheppard, R. Terrell and G. Henkelman, J. Chem. Phys., 2008, 128, 134106.
76. S. Sadan, I.-H. Svenum, S. Ø. Hanslin and J. Akola, Phys. Chem. Chem. Phys., 2024, 26, 25957–25968.
77. M. Hakala and K. Laasonen, Phys. Chem. Chem. Phys., 2018, 20, 13785–13791.
78. S. S. Yadavalli, G. Jones and M. Stamatakis, Phys. Chem. Chem. Phys., 2021, 23, 15601–15612.
79. L. Zhu, C. Liu, X. Wen, Y.-W. Li and H. Jiao, Catal. Sci. Technol., 2019, 9, 199–212.
80. K. Tonigold and A. Groß, J. Comput. Chem., 2012, 33, 695–701.
81. Q. Li and X. Hu, Phys. Rev. B:Condens. Matter Mater. Phys., 2006, 74, 035414.
82. R. B. Wexler, J. M. P. Martirez and A. M. Rappe, Chem. Mater., 2016, 28, 5365–5372.
83. J. He, Á. Morales-García, O. Bludsky and P. Nachtigall, CrystEngComm, 2016, 18, 3808–3818.
84. H. Jónsson, G. Mills and K. W. Jacobsen, Classical and quantum dynamics in condensed phase simulations, World Scientific, 1998, pp. 385–404.
85. G. Henkelman and H. Jónsson, J. Chem. Phys., 2000, 113, 9978–9985.
86. G. Henkelman, B. P. Uberuaga and H. Jónsson, J. Chem. Phys., 2000, 113, 9901–9904.
87. P. Xiao, D. Sheppard, J. Rogal and G. Henkelman, J. Chem. Phys., 2014, 140, 174104.
88. J. Kästner and P. Sherwood, J. Chem. Phys., 2008, 128, 014106.
89. A. Heyden, A. T. Bell and F. J. Keil, J. Chem. Phys., 2005, 123, 224101.
90. H. Group, VTST Tools, https://theory.cm.utexas.edu/vtsttools/scripts.html, Accessed: 16 September 2025.
91. K. Mathew, V. S. C. Kolluru and R. G. Hennig, VASPsol: Implicit solvation and electrolyte model for density-functional theory, https://github.com/henniggroup/VASPsol, 2018, https://github.com/henniggroup/VASPsol.
92. K. Mathew, R. Sundararaman, K. Letchworth-Weaver, T. Arias and R. G. Hennig, J. Chem. Phys., 2014, 140, 084106.
93. K. Mathew, V. Kolluru, S. Mula, S. N. Steinmann and R. G. Hennig, J. Chem. Phys., 2019, 151, 234101.
94. S. Petrosyan, J.-F. Briere, D. Roundy and T. Arias, Phys. Rev. B:Condens. Matter Mater. Phys., 2007, 75, 205105.
95. G. Henkelman and H. Jónsson, J. Chem. Phys., 1999, 111, 7010–7022.
96. K. Chan and J. K. Nørskov, J. Phys. Chem. Lett., 2015, 6, 2663–2668.
97. R. E. Warburton, A. V. Soudackov and S. Hammes-Schiffer, Chem. Rev., 2022, 122, 10599–10650.
98. S. Ghosh, A. V. Soudackov and S. Hammes-Schiffer, J. Chem. Theory Comput., 2016, 12, 2917–2925.
99. M. Kessinger, A. V. Soudackov, J. Schneider, R. E. Bangle, S. Hammes-Schiffer and G. J. Meyer, J. Am. Chem. Soc., 2022, 144, 20514–20524.
100. P. Hutchison, R. E. Warburton, A. V. Soudackov and S. Hammes-Schiffer, J. Phys. Chem. C, 2021, 125, 21891–21901.
101. M.-H. Ho, R. Rousseau, J. A. Roberts, E. S. Wiedner, M. Dupuis, D. L. DuBois, R. M. Bullock and S. Raugei, ACS Catal., 2015, 5, 5436–5452.
102. I. Ledezma-Yanez, W. D. Z. Wallace, P. Sebastián-Pascual, V. Climent, J. M. Feliu and M. Koper, Nat. Energy, 2017, 2, 1–7.
103. D. Guo, Y. Nakagawa, H. Ariga, S. Suzuki, K. Kinoshita, T. Miyamoto, S. Takakusagi, K. Asakura, S. Otani and S. T. Oyama, Surf. Sci., 2010, 604, 1347–1352.
104. A. B. Hernandez, H. Ariga, S. Takakusagi, K. Kinoshita, S. Suzuki, S. Otani, S. T. Oyama and K. Asakura, Chem. Phys. Lett., 2011, 513, 48–52.
105. P. Liu, J. A. Rodriguez, T. Asakura, J. Gomes and K. Nakamura, J. Phys. Chem. B, 2005, 109, 4575–4583.
106. S. D. Ghadge, O. I. Velikokhatnyi, M. K. Datta, P. M. Shanthi and P. N. Kumta, Catal. Sci. Technol., 2021, 11, 861–873.
107. S. Sadan, HER on Ni2P catalyst voltage calculation data, Zenodo, 2026 DOI:10.5281/zenodo.18991779.

---