An efficient and cyclic hydrogen evolution reaction mechanism on [Ni(P^H₂N^H₂)₂]²⁺ catalysts: a theoretical and multiscale simulation study

Abstract

In this paper we report a theoretical and a multiscale simulation study of the hydrogen evolution reaction (HER) on the [Ni(P^H₂N^H₂)₂]²⁺ catalyst in acidic media (2H⁺ + 2e⁻ → H₂). First, at the DFT calculations level, a cyclic pathway for the HER is proposed highlighting the shuttling of electrons with protons on the conformationally flexible catalyst. The theoretical calculation gives a better understanding of the efficient cyclic pathway of [Ni(P^H₂N^H₂)₂]²⁺, and the effect of solvent on the mechanism has been discussed. The σ-donating and π-accepting nature of H₂–Ni bond has been identified in the H₂ complex. The oxidation state of the Ni centre and geometrical changes of the catalyst in the reaction coordinate are also identified. Then a mean-field kinetic model incorporating the calculated DFT data has been developed. This model allows us to simulate the behaviour of these catalysts in electrochemical conditions representative of polymer electrolyte membrane water electrolyzers operation. Calculated results include experimental observables such as polarization curves showing good agreement with available experimental data. Competitive phenomena between the different electrochemical mechanisms, the protons and H₂ transport, and their relative impact on the overall cell performance are particularly discussed.

---

Introduction

The production of molecular hydrogen (H₂) is a main step towards the widespread fuel cell development as an electric power generator. Thus, the efficient generation, storage, and oxidative conversion of H₂ into electric power via fuel cells have been the theme of intense studies.^1,2 In recent years platinum based materials have attracted much research attention as proton reductive catalysts in polymer electrolyte membrane water electrolyzer (PEMWE) cathodes.³ Some serious challenges still need to be overcome before the development and commercialisation of these technologies, such as their high catalyst cost and durability.⁴

In biology, various iron hydrogenase are known to be involved in reductive generation and oxidative uptake of molecular hydrogen.⁵ Fe-only hydrogenase,⁶ Ni–Fe hydrogenase⁷ and Fe–Fe hydrogenase⁸ are well-known examples of enzymes catalyzing the hydrogen oxidation reaction (HOR) and the hydrogen evolution reaction (HER). Ni–Fe hydrogenase has been demonstrated to produce molecular hydrogen efficiently and this has received substantial attention in recent years.^9–12 The bimetallic active site and outer-coordination sphere of this enzyme is believed to control the efficient formation or cleavage of the H–H bond.¹³

Numerous functional model complexes have been identified and reported experimentally to mimic the HER functionality of Ni–Fe hydrogenases.^14,15 Le Goff et al., have reported a new type of model nickel bisdiphosphine ([Ni(P₂^HN₂^H)₂]²⁺) to mimic the HER function of Ni–Fe hydrogenase enzyme.¹⁶ This system exhibits electrochemical activity towards the production of molecular hydrogen from protons and electrons (2H⁺ + 2e⁻ → H₂). The general reversible cycle mechanism of H₂ production for the Ni(P₂^HN₂^H)₂]²⁺ catalyst is presented in Fig. 1.

Fig. 1 Schematic representation of reversible hydrogen evolution reaction of Ni(P₂HN₂H)₂]²⁺.

In a previous publication, we have studied the H₂ production mechanism from the proton reduction the [Ni(P₂^HN₂^H)₂]²⁺ catalyst, concluding that the [H₂–Ni(P₂^HN₂^H)₂]²⁺ complex is formed via a symmetrical dissociation of N–H bond from molecule A.¹⁷ One year later, more insight has been gained by Dupuis et al., who have reported that such symmetric TS was not related to the reaction mechanism. They concluded that the H₂ oxidation or proton reduction with the [Ni(P₂^HN₂^H)₂]²⁺ family of electrocatalysts proceeds through a heterolytic cleavage of N–H bond of A molecule and formation of [H₂–Ni(P₂^HN₂^H)₂]²⁺ complex and a metal hydride proton species that can evolve further into diproton species.¹⁸

In this context, two issues are addressed in this paper. Firstly a complete cyclic pathway for molecular HER of [Ni(P₂^HN₂^H)₂]²⁺ catalyst is proposed and studied in more details using DFT. In particular the electronic structure and solvent effects along the HER are investigated. Secondly, in order to validate the proposed mechanism of HER and the electrochemical behaviour of [Ni(P₂^HN₂^H)₂]²⁺ in acidic media, a mean-field kinetic model, following the multiscale approach previously developed by Franco¹⁹ and others^20–22 has been developed using DFT-calculated thermodynamic parameters and experimentally fitted parameters. The multiscale simulation results are finally compared with the experimental data.

Computational details

All structures were optimized at the density functional theory (DFT) level with the PBE^23,24 density functional and a mixed basis set which was made up of the Stuttgart Dresden relativistic ECP basis set (SDD)²⁵ for Ni and the 6-31G** basis set for all non-metal atoms. The transition states, intermediates and products are calculated at B3LYP/SDD level using Gaussian 03²⁶ and stationary points have been characterized by computing vibrational frequencies. Transition states have been further confirmed by harmonic frequency analysis and by IRC and observation of single imaginary frequency. NBO analysis has been carried out for selected transition states geometries at the same level of theory.^27–30 Solvent effects (water and acetonitrile, chosen because they are experimentally considered) were modelled using the polarized-continuum-model^31–33 (PCM), on the B3LYP/SDD optimized structures. Ni²⁺ (3d⁸) complexes were always considered as low spin species as expected in square planar symmetry.

Results and discussion

Cyclic reaction pathways and energetics of catalytic H₂ evolution mechanism of [Ni(P₂^HN₂^H)₂]²⁺

The cyclic pathway for reversible H₂ evolution of [Ni(P₂^HN₂^H)₂]²⁺ complex is presented in Scheme 1. The corresponding calculated relative free energy profile is shown in Fig. 2. All cartesian coordinates of optimized geometries of reactants, transition states, intermediates and products are given in the ESI.†

Fig. 2 Relative free energy profiles for HER of Ni(P₂^HN₂^H)₂]²⁺ computed at B3LYP/SDD level in gas phase, and using a PCM model for solvation using H₂O (red colour dotted line) and CH₃CN (green colour dotted line).

Initially the optimized geometry of [Ni(P₂^HN₂^H)₂]²⁺ catalyst is square planar with Ni²⁺ oxidation state. This complex readily accepts two electrons and two protons from the solvent and electrode. Unfortunately we cannot have any estimation of the energy barrier of such a step as it involves molecules with different numbers of atoms. So that in further development of the kinetic model, we made the assumption that the activation barrier is zero and this step is just governed by the concentration of protons and the interaction with the electron providing interface of the electrodes. The protons coupled to electrons transfer to [Ni(P₂^HN₂^H)₂]²⁺ from acid source leads to the diprotonated species (A) and the optimized geometry of A reveals a tetrahedral coordination with Ni⁰ oxidation state. The calculated occupancy (Table 1) of metal orbitals (4s = 1.50 and 3d = 7.3) is consistent with a zero oxidation state of nickel (3d⁸4s² atomic configuration). This step is not shown in Fig. 2.

Table 1 The occupancy of 3d and 4s orbital of A, B and C species

Reactive species	Occupancy
	3d	4s
A	7.3	1.50
B	7.3	0.82
C	7.4	0.10

---

First intramolecular proton transfer step (A → B)

The diprotonated A molecule leads to a nickel hydride intermediate (B) through intramolecular proton transfer transition state (TS_AB) with an activation energy of 9.2 kcal mol⁻¹ in gas phase. The intramolecular proton transfer transition state is further confirmed by Intrinsic Reaction Coordinate (IRC) calculation giving a single imaginary frequency. The evolution from A to the TS_AB of molecule A is accompanied by a change from a tetrahedral to a distorted square pyramidal symmetry around Ni and an increase of the oxidation state of Ni from 0 to 1+. The occupation of the 4s orbital is decreased from 1.50 to 0.82 during this step, in agreement with a change in nickel oxidation state from 0 to 1.

Second intramolecular proton transfer step (B → C)

The nickel hydride (B) undergoes another one intramolecular proton transfer (TS_BC) reaction leading to the H₂–[Ni(P₂^HN₂^H)₂]²⁺ complex (C) with an activation energy of 16.1 kcal mol⁻¹ in gas phase. The formation of H₂ complex (C) through TS_BC is the rate determining step due to the requirement of higher activation energy. The electron occupancy of 4s and 3d orbitals (resp. 0.10 and 7.4) of Ni in C complex is lower than the occupancy of 4s orbital (0.82) in B molecule and corroborates the change of nickel oxidation state from 1+ to 2+.

Hydrogen evolution step (C → D)

The H₂ complex (C) undergoes H₂ evolution reaction through symmetrical breaking of η²_Ni–H2 bond (TS_CD) and gives a square planar geometry with Ni²⁺ oxidation state D molecule as product. The required activation energy for breaking the η²_Ni–H2 coordination bond to Ni has been calculated as 14.5 kcal mol⁻¹ in gas phase. In addition, a rotational transition state (TS_CC) is also observed from H₂–[Ni(P₂^HN₂^H)₂]²⁺ complex at 13.4 kcal mol⁻¹ in gas phase. But the IRC calculation of the TS_CC leads to the same structure as reactant C with the same energy and this step has not been included in the relative free energy profile for clarity purpose.

Whether the catalyst is involved in the cyclic pathway or interacts with the solvent?

After the formation of molecule D, two types of reaction have been considered here (Scheme 1). One is the interaction of D molecule with acetonitrile¹⁶ solvent leading to the solvated product (E). Experimentally, the solvent is acetonitrile and we have chosen the same as ligand in our computation. The molecule D may interact with acetonitrile solvent through TS_DE to give an acetonitrile bonded product (E) which is thermodynamically favoured. The interaction of solvent molecule requires very high activation energy (16.2 kcal mol⁻¹). Secondly, by a new 2e⁻ and 2H⁺ transfer, the molecule D leads to the kinetically more favourable initial product (A).

Scheme 1 Schematic representation of various steps involving reversible hydrogen evolution reaction of Ni(P₂^HN₂^H)₂]²⁺.

Bulk solvents effect on activation energies of HER

In order to check the influence of solvent on the reaction pathway, the polarizable continuum model has been applied to the calculations of the free energies of the reaction pathway. Two solvents have been considered, water and acetonitrile, and the calculated solvation free energies are given in the profile in Fig. 2. The reaction profile shows increases in activation and free energies of respectively the transition states and intermediates with respect to the gas phase. This solvent destabilization was expected due to the screening of electrostatic interaction in a medium with a high dielectric constant. The destabilization effect of H₂O is slightly higher compared to CH₃CN towards the reaction process. We can conclude that the effect of solvent is an upward shift of all energies, but does not alter the reaction pathway and geometries.

Changes of geometry and oxidation states by intra molecular electrons migration

The central Ni atom undergoes several geometry and charge changes along the reaction coordinates. In order to correlate the geometrical changes and variation of Ni charge, the NBO analysis and d-orbitals occupation and energy estimation have been done for all transition states and intermediate species.

Firstly we have examined the relation between Ni^0/2+ oxidation state and optimized geometries for Ni²⁺ and Ni⁰ respectively in the species [Ni(P₂^HN₂^H)₂]²⁺ (D) and [Ni⁰(P₂^HN₂^HH)₂]²⁺(A), the latter representing the diprotonated complex. The optimization profiles of both cases are given in ESI (S1†). As expected, the Ni²⁺ state of nickel ([Ni(P₂^HN₂^H)₂]²⁺) always prefers a tetrahedral geometry while Ni⁰ ([Ni(P₂^HN₂^HH)₂]²⁺) always prefers a square planar geometry. This strong correlation between the coordination geometry and oxidation state is thus a very good tool in the determination of the oxidation state of nickel atom in the intermediate species.

The energy splittings of Ni 3d orbitals in each step of the reaction have been computed and presented in Fig. 3. The optimized geometry of A gives the following energy levels of metal d-orbitals: d_xy (−10.3 eV) ≈ d_yz (−10.5 eV) ≈ d_zx (−10.5 eV) > d_x2–y2 (−11.8 eV) ≈ d_z2 (−11.7 eV), in agreement with a tetrahedral symmetry around nickel.

Fig. 3 The energy variation profile of metal orbitals along with geometrical changes in the reaction coordinates.

In the TS_AB transition state, the energy splitting of d orbitals is d_z2 (−9.5 eV) > d_xy (−10.5 eV) ≈ d_x2–y2 (−10.5 eV) > d_xz (−11.2 eV) ≈ d_yz (−11.3 eV) consistent with the trigonal bipyramidal geometry. The high lying metal orbital is d_z2 due to the presence of ligands on z axis in the movement of A to TS_AB. In order to quantify the oxidation state changes of Ni, the NBO analysis has been carried out on the TS_AB optimized geometry. The charge transfer probability of d(Ni_dxz) → σ* (N–H) interaction is 56.6 kcal mol⁻¹ in TS_AB. The transfer of electrons from bonding level d(Ni_dxz) to anti bonding σ* (N–H) level decreases the N–H bond order and increases the oxidation state of Ni from 0 to 1+. The occupancy of d_xz orbital in TS_AB is comparatively lower than that of d_xz orbital in A.

Then during the change from TS_AB to the nickel hydride B molecule, the high-lying d_z2 orbital is stabilized and the low-lying d_x2–y2 orbital is destabilized. The 3d orbital splitting of nickel in B is d_x2–y2 (−9.8 eV) > d_z2 (−10.1 eV) > d_xy (−11.0 eV) > d_xz (−12.7 eV) ≈ d_yz (−12.6 eV) in agreement with the square pyramidal geometry. The NBO analysis shows that the occupation of d_x2–y2 orbital is less than two. This observation confirms the change of Ni oxidation state from 0 to 1+ during the TS_AB and this result exactly matches recently reported electrochemical studies on similar Ni model complexes.³⁴

In order to form the H₂–[Ni(P₂^HN₂^H)₂] C complex, the molecule B undergoes another intramolecular proton movement through the transition state (TS_BC), by transfer of electron density from d_xz orbital to σ* of N–H bond. The energy level splitting pattern of d-orbital of nickel in TS_BC is d_x2–y2 (−11.0 eV) > d_z2 (−11.5 eV) > d_xy (−12.8 eV) > d_xz (−13.1 eV) ≈ d_yz (−13.2 eV), consistent with the square pyramidal geometry. The transferred proton from nitrogen is simultaneously connected with H atom and Ni atom (Scheme 1) described after and finally TS_BC leads to the (η²-H₂)-[Ni(P₂^HN₂^H)₂]²⁺ (C) complex. The second order perturbation energy calculated for d(Ni_dxz2) → σ*(N–H) interaction in TS_BC, is 34.4 kcal mol⁻¹. Moreover, the occupancy of d_z2 of nickel in TS_BC is comparatively lower than the occupancy of d_z2 orbital in B molecule. These observations clearly indicate that two electrons are transferred from Ni atom to two σ*(N–H) orbitals for each intramolecular transition states (TS_AB & TS_BC) and these two electrons are “consumed” for H–H bond formation. The product of TS_BC transition state is the H₂–[Ni(P₂^HN₂^H)₂]²⁺ (C) complex with the H₂ molecule connected to Ni atom by an η²-H₂ coordination fashion. This η²-H₂ coordination mode yields a distorted square pyramidal geometry. The calculated oxidation state changes of nickel during this intramolecular electron transfer in the transition state nicely agrees with recent report of a similar model complex.³⁴

The energy level splitting pattern of 3d orbitals of nickel in C is d_z2 (−11.4 eV) > d_x2−y2 (−11.8 eV) > d_xy (−12.7 eV) > d_xz (−13.4 eV) ≈ d_yz (−13.5 eV) consistent with a distorted square pyramidal geometry. A 3 center-2 electron bonding scheme is observed in this (η²-H₂)-[Ni(P₂^HN₂^H)₂]²⁺ complex (Fig. 4). It is stabilized by some charge transfers that have been revealed by the NBO analysis. A σ-donation from H–H bonding orbital electron density to d_z2 orbital of nickel is observed, as well as π back-donation from d_xz (Ni) to σ*(H–H). The optimized geometry of (η²-H₂)-[Ni(P₂^HN₂^H)₂]²⁺ shows a 0.8 Å length for H–H bond and this observation nicely agrees with experimental report for true H₂–M complex.³⁵ The perturbation energy of σ donating σ_(H–H) → d_z2 (Ni) interaction is 23.9 kcal mol⁻¹ and π back bonding d_xz (Ni) → σ*_(H–H) interaction energy is 20.5 kcal mol⁻¹. The perturbation energy of σ donor interaction is relatively higher compared to π back bonding.

Fig. 4 Pictorial views of π-acceptor (1) and σ-donor electron transfer of H₂–Ni(P₂^HN₂^H)₂]²⁺ complex.

The highly distorted square pyramidal geometry of (η²-H₂)–[Ni(P₂^HN₂^H)₂]²⁺ complex undergoes H₂ evolution through TS_CD. The splitting of d-orbitals energies in TS_CD is very similar to the d-orbital splitting of C, but the individual energy levels increase in the TS_CD. The splitting pattern of orbital energy level in TS_CD is d_z2 (−9.1 eV) > d_x2–y2 (−10.5 eV) > d_xy (−11.7 eV) > d_xz (−12.9 eV) ≈ d_yz (−12.8 eV) consistently with the square pyramidal geometry.

In the generation of the H₂ species bonded in η² fashion to Ni (step C to D), it is interesting to investigate the structural parameters associated to H₂, in the complexes involved in this process. The Ni–H bond distance (1.93 Å) in (η²-H₂)-[Ni(P₂^HN₂^H)₂]²⁺ is relatively lower compared to the Ni–H bond distance (2.10 Å) of the TS_CD species. Conversely, the H–H distance (0.79 Å) of η²-H₂–[Ni(P₂^HN₂^H)₂]²⁺ is relatively higher compared to H–H distance (0.77 Å) of TS_CD. These observations clearly point out the elongation of Ni–H bonds and shortening of H–H bond during the moving of C to TS_CD, that will lead further to the release of the H₂ molecule.

The H₂ evolution step (TS_CD) in the pathway of HER leads to the square planar D molecule with Ni²⁺ oxidation state. The 2+ oxidation state of nickel is confirmed by occupancy estimation (4s = 0.10 and 3d = 7.3) of 4s orbital in D molecule by NBO analysis. The d-orbital splitting pattern of D molecule is d_x2–y2 (−9.4 eV) > d_xy (−10.3 eV) > d_xy (−11.7 eV) > d_xz (−12.2 eV) ≈ d_yz (−12.3 eV) and this arrangement is consistent with the square planar geometry of molecule D. The geometrical features of TS_CC, TS_DE and E species in the pathway all correspond to square pyramidal geometry with the same d-orbital splitting pattern.

To summarize this section, the NBO and individual d-orbitals energy level estimation offer nice insights into the geometrical and oxidation state changes of the [Ni(P₂^HN₂^H)₂]²⁺ catalyst along with reaction coordinates. The initial oxidation state 2+ of nickel is converted to 0 by accepting two electrons and two protons in the initial step. The TS_AB and TS_BC result in a 1+ oxidation state of nickel by the transfer of two electrons from high-lying d-orbital to σ*(N–H) orbital of deprotonating nitrogen (Scheme 1). The σ-donation and π-back donation of H₂–Ni bond has been identified in the H₂ complex. These analyses give some clues on the origin of the formation of the H₂ bond. Our calculated pathway is consistent with part of Dupuis et al., report¹⁸ unravelling the influences of geometrical, electronic and solvent effects and finally it completes the cyclic mechanism.

Mean field kinetic model

The multiscale simulation approach developed by Franco et al. allows describing the detailed physicochemical processes at multiple scales in electrochemical devices for energy conversion and storage, such as PEMFC, electrolyzers and batteries.^19,36–39 Such a model is a multiscale one in the sense that it is made of a set of interconnected sub-models describing the phenomena occurring at different scales. The model results from the coupling of different geometrical scales: e.g. reactants and charge 1-D mesoscale description across the membrane-electrodes assembly (MEA) and the interfacial nanoscale mechanisms at the vicinity of the electrochemically active materials including both elementary kinetics (with ab initio calculated activation barriers) and electrochemical double layer effects. The description of the mechanisms in this approach remains macroscopic in the sense that it is based on irreversible thermodynamic concepts as they are extensively used in chemical engineering: use of conservation laws coupled to closure equations (flux expressions, chemical rate models, and thermodynamic models). The model is implemented in a Matlab/Simulink environment coupled with numerical modules programmed in C language.³⁷

In this paper we have adapted this model for describing the HER on the [Ni(P₂^HN₂^H)₂]²⁺ catalysts supported on Multi-Walled Carbon Nanotubes (MWCNTs).

Global model structure

The real PEMWE cathode is represented as being constituted by three coupled spatial scales (Fig. 5): one describing the electrochemical double layer behaviour under non equilibrium (nanoscale), one describing the transport of H₂ and protons across the electrolyte film at the vicinity of the catalyst (microscale) and one describing the overall electrode (macroscale). The electrochemical double layer model is further divided in two layers: the compact and the diffuse layers.³⁶ The bio-inspired catalyst material is assumed to be linked in turn to the MWCNTs supported on an Indium Tin Oxide (ITO) substrate.¹⁶

Fig. 5 Schematic representation of developed fuel cell model and nanoscale representation of electrical double layer including bioinspired nanomaterial coated compact layer.

The constructed numerical code allows the analysis of the stationary and dynamic behaviour of the different state variables (for example, C_H⁺(t), σ(t), etc.), the electrode electronic potentials, ψ(t) in response to a current demand, I(t), at a given temperature, T, and pressure P_cathode.

Microscale model

The microscale model describes H₂ transport across the electrolyte film thickness. The H₂ product is assumed to diffuse mainly perpendicularly to the catalytic layer. The electrons are transported to the catalyst molecules from the external electrical circuit.

The flux for the electrically neutral species H₂ is assumed to be given by a Fick's law of diffusion and the mass balance leads to

where D_H2 is the diffusion coefficient in the electrolyte assumed to be constant.

Nanoscale model (electrochemical double layer)

The pictorial representation of the electrochemical double layer has been given in Fig. 6. The catalyst molecules are imagined to be aligned and an imaginary reaction plane is placed at x = L. Adsorption of reaction intermediates, water and products (forming the compact layer) is considered to take place on this imaginary plane. The detailed description of this model is given below.

Fig. 6 Schematic representation of nanoscale model for electrochemical cathode in electrode of fuel cell.

Diffuse layer

The diffuse layer model describes the transport by diffusion and migration of the reaction species coupled with the electric field generated by the resulting charge distribution (Fig. 6). The H₂, H⁺ and H₂O species are considered as punctual. The solvation and convection by water are not considered here.

In a similar way that in the H₂ case, combining the flux related to the Fick's diffusional force to the electrical force with the mass balance, the following equations have been derived representing the H⁺ transport:

where ϕ(x, t) is the electric potential in the electrolyte, which is calculated from the Poisson's equation.¹⁹

The scalar potential ϕ(x = 0, t) equals 0 at x = L, the boundary conditions for the mass balance equations for the neutral and charged species are

J_H2 = (L, t) = (v_TAF + v_HEY) (3)

(positive because H₂ is produced)(negative because protons are consumed).

The kinetic rates in eqn (4) and (5) are provided below.

Compact layer

Within the elementary kinetic model, the following global reaction is considered:

[Ni(P^H₂N^H₂)₂]²⁺ + CH₃CN + 2e⁻ + 2H⁺ → [Ni(P^H₂N^H₂)₂(CH₃CN)]²⁺ + H₂ (5)

The kinetic parameters (associated with each elementary kinetic step) are assumed to be given by

where ΔG^‡ is the activation barrier at room temperature. These activation barriers are the ones reported in Fig. 2. The k_B and h are the Boltzmann and Planck's constants respectively.

f(η_F(t)) is a function of the electrostatic potential jump between the catalyst surface and the electrolyte, through the compact layer and accounts for the effect of the interfacial electric field on the electrochemical reactions.³⁶

According to the kinetic model, the molecular hydrogen (H₂) generation mechanism is fully depending on the Ni and N active sites. The coverage balance equation of Ni active sites can be written as

where θ^Ni_H and are the coverage of reaction intermediates, is the coverage by the CH₃CN, and and are the water coverage following up (pointing out) and down (towards) orientation to the Ni site.

The coverage balance equation for the N sites in the catalyst is

where θ^N is the free sites coverage and θ^N_H is the reaction intermediate coverage.

According to the proposed mechanism of HER of [Ni(P₂^HN₂^H)₂]²⁺, the whole reactions have been divided into six different elementary kinetic steps (refer to Scheme 1) given below: kinetics expression of initial electron coupled protonation of catalyst.

v₁ = k₁(θ^N)²(a_H⁺)²e^−2αfn − k₋₁(θ^N_H)²e^2(1−α)fn (9)

kinetics of first intramolecular proton transfer (A → B)

v₂ = k₂(θ^N_H)²(θ^Ni) − k₋₂(θ^Ni_H)(θ^N_H)(θ^N) (10)

kinetics of second intramolecular proton transfer (B → C)kinetics of rotational transition state (C → C)kinetics of H₂ evolution step (C → D)kinetics of solvent absorption of catalyst (D → E)

The surface potential η^F is defined as the electrostatic potential difference across the compact layer.³⁶

The coverage of all nickel and nitrogen reaction intermediates are obtained by solving the following set of equations.

The charge density σ(t) has been calculated from the current density conservation at the catalyst/electrolyte interface according to:

where J(r, t), is the local electronic current density traversing the catalysts layer. For our bio-inspired molecules, we demonstrated that the surface potential is given by:where σ(t) is the surface electronic charge density and X is a transcendental function of the coverage fractions in intermediate reaction species and the charge density as have been reported by Franco et al.^35–38

In equation (16), the parameter a_N is given by:

where ΔG_equ is the water adsorption energy given by:

ΔG_equ = E_Catalyst+H2O − E_Catalystfree − E_H2O (18)

Simulation results

After deriving all the elementary steps involved in the H₂ production mechanism towards the Pt-free catalyst, the observable properties of current versus electrode potential have been calculated using the thermodynamical parameters of each elementary kinetic steps which have been calculated in DFT. The comparison of simulated and experimental current density versus potential curves has been made and the result is presented in Fig. 7. The simulated current density versus potential change curve is in reasonable agreement with experimentally observed data.¹⁶

Fig. 7 Evolution of current density as a function of the potential for the HER by Ni(P₂^HN₂^H)₂]²⁺ catalyst. The solid line indicates experimental value and dotted line indicates the simulated curve.

The coverage changes of reaction intermediates have been calculated as a function of potential variation towards the HER and presented in Fig. 8–10. The changes of nitrogen and nickel coverage with respect to potential changes have been illustrated in Fig. 9 and 10. The coverage of nitrogen intermediate (θ^NiH) is decreased along with potential increase (i.e. when potential becomes more positive) while coverage of nitrogen free sites (θ^N) increases as the potential increases. A very slight variation is observed for the nickel intermediates coverage telling us that the reactions on the nitrogen sites are governing the overall electrochemical activity of the catalysts towards the HER. These two species are involved in the initial diprotonation step of [Ni(P₂^HN₂^H)₂]²⁺ catalyst and the simulation results reveal that initial diprotonation step is responsible for electrochemical HER. This observation is consistent with experimental knowledge.¹⁶

Fig. 8 Changes of coverage of reaction intermediates along with variation of potential in the hydrogen evolution reaction of Ni(P₂^HN₂^H)₂]²⁺ catalyst.

Fig. 9 Changes of coverage of nitrogen and nitrogen intermediate along with variation of potential in the hydrogen evolution reaction of Ni(P₂^HN₂^H)₂]²⁺ catalyst.

Fig. 10 Influence of water on the electrochemically accessible Ni free sites on the electrode surface during hydrogen evolution reaction catalysed by Ni(P₂^HN₂^H)₂]²⁺.

The detail of the influence of water molecules on the electrochemically accessible area of Ni free site of the electrode during the HER catalysed by [Ni(P₂^HN₂^H)₂]²⁺ is presented in Fig. 10. The coverage change with potential is small but the proportion of up and down dipoles is controlled by the charge density of the catalyst.

Conclusions

DFT calculations have been used to investigate the HER catalysed by [Ni(P₂^HN₂^H)₂]²⁺. The hydrogen evolution mechanism of [Ni(P₂^HN₂^H)₂]²⁺ catalyst has been considered as involving six different steps in a cyclic pathway and all the steps are connected through transition states. The geometrical changes of the catalyst along with reaction pathway have been identified. The oxidation state of nickel is changed from 2+ to 0 during the course proton absorption. Furthermore the oxidation state of nickel atom changed from 0 to 1+ and 1+ to 2+ during the H₂ molecule evolution reactions and these changes are nicely brought out from the NBO analysis. The observed oxidation state and geometrical changes of catalyst are consistent with recently reported experimental results. The σ-donation and π-back donation contributions in H₂–Ni bond have been identified in the H₂ complex. It gives important information about the 2e⁻ transfer necessary for the H₂ bond.

Understanding of elementary kinetic steps of mechanisms of electrochemical reactions is crucial for the identification of those steps that are the major contributors to overpotential losses. Furthermore the multiscale simulation approach of electrochemical power generators, developed by Franco et al., has been used here to simulate the features of the HER on these catalysts, using DFT computed thermodynamical data. The calculated polarization curve shows reasonable agreement with available experimental data. These simulations allowed us to find the reaction intermediate concentration and the effect of electrochemically accessible protonating sites. The variations of nitrogen free site and protonating nitrogen site are linearly dependant on the variation of potential in the simulated system. The cyclic mechanistic pathway allows to understand the efficient hydrogen evolution nature of [Ni(P₂^HN₂^H)₂]²⁺ catalyst at the molecular level. Apart from the fundamental interests in understanding these systems, this modelling approach would allow also to consider amendment of the catalyst in order to improve its performance.

Acknowledgements

The authors thank the CEA Nanosciences program for financial support. SD dedicated to Professor P. Venuvanalingam on his 60th Birthday. SD wishes to thank the European Community for a EUROTALENTS grant, from the COFUND Marie Curie actions in FP7.

References

1. E. S. Wiedner, J. Y. Yang, W. G. Dougherty, W. S. Kassel, R. M. Bullock, M. R. DuBois and D. L. DuBois, Organometallics, 2010, 29, 5390–5401.
2. S. Horvath, L. E. Fernandez, A. V. Soudackov and S. Hammes-Schiffer, Proc. Natl. Acad. Sci. U. S. A., 2012, 15663–15668.
3. F. Barbir, Sol. Energy, 2005, 78, 661–669.
4. A. A. Franco, PEMFC degradation modeling and analysis, in Polymer electrolyte membrane and direct methanol fuel cell technology (PEMFCs and DMFCs) Volume 1: Fundamentals and performance, ed. C. Hartnig and C. Roth, Woodhead, Cambridge, UK, 2012.
5. A. Ciaccafava, P. Infossi, M. Ilbert, M. Guiral, S. Lecomte, M. T. Giudici-Orticoni and E. Lojou, Angew. Chem., Int. Ed., 2012, 51, 953–968.
6. J. W. Peters, W. N. Lanzilotta, B. J. Lemon and L. C. Seefeldt, Science, 1998, 282, 1853–1858.
7. R. P. Happe, W. Roseboom, A. J. Plerik, S. P. J. Albracht and K. A. Bagley, Nature, 1997, 385, 126.
8. J. W. Tye, M. P. Hall and M. Y. Darensbourg, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 16911–16912.
9. Hydrogen as a Fuel: Learning from Nature, ed. R. Cammack, M. Frey and R. Robson, Taylor & Francis, London, 2001.
10. M. A. Alonso-Lomillo, O. Rüdiger, A. Maroto-Valiente, M. Velez, I. R. guez-Ramos, F. J. Muňoz, V. M. Fernández and A. L. De Lacey, Nano Lett., 2007, 7, 1603–1608.
11. M. Frey, ChemBioChem, 2002, 3, 153–160.
12. Y. Nicolet, A. L. de Lacey, X. Vernede, V. M. Fernandez, E. C. Hatchikian and J. C. Fontecilla-Camps, J. Am. Chem. Soc., 2001, 123, 1596–1601.
13. A. Jain, S. Lense, J. C. Linehan, S. Raugei, H. Cho, D. L. DuBois and W. J. Shaw, Inorg. Chem., 2011, 50, 4073–4085.
14. C. Tard and C. J. Pickett, Chem. Rev., 2009, 109, 2245–2274.
15. A. Fihri, V. Artero, M. Razavet, C. Baffert, W. Leibl and M. Fontecave, Angew. Chem., 2008, 120, 574–577.
16. A. Le Goff, V. Artero, B. Jousselme, P. D. Tran, N. Guillet, R. Métayé, A. Fihri, S. Palacin and M. Fontecave, Science, 2009, 326, 1384–1387.
17. A. Kachmar, V. Vetere, P. Maldivi and A. A. Franco, J. Phys. Chem. A, 2010, 114, 11861–11867.
18. M. Dupuis, S. Chen, S. Raugei, D. L. DuBois and R. M. Bullock, J. Phys. Chem. A, 2011, 115, 4861–4865.
19. A. A. Franco, P. Schott, C. Jallut and B. Maschke, Fuel Cells, 2007, 07, 99–107.
20. Y. Yang and S. Holdcroft, Fuel Cells, 2005, 5, 171–186.
21. M. A. Hickner and B. S. Pivovar, Fuel Cells, 2005, 5, 213–229.
22. H. E. Zhang and Z. T. Zhou, Prog. Chem., 2008, 20, 602.
23. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
24. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1997, 78, 1396.
25. D. Andrae, U. Häussermann, M. Dolg, H. Stoll and H. Preuss, Theor. Chem. Acc., 1990, 77, 123–141.
26. M. J. Frisch, et al., Gaussian 03, revision E.01, Gaussian, Inc., Wallingford, CT, 2004.
27. J. P. Foster and F. Weinhold, J. Am. Chem. Soc., 1980, 102, 7211–7218.
28. A. E. Reed and F. Weinhold, J. Chem. Phys., 1983, 78, 4066–4073.
29. A. E. Reed, R. B. Weinstock and F. Weinhold, J. Chem. Phys., 1985, 83, 735–746.
30. A. E. Reed and F. Weinhold, J. Chem. Phys., 1985, 83, 1736–1740.
31. M. C. F. Wander and A. E. Clark, Inorg. Chem., 2008, 47, 8233–8241.
32. B. Mennucci and J. Tomasi, J. Chem. Phys., 1997, 107, 3032–3041.
33. B. Mennucci, E. Cances and J. Tomasi, J. Phys. Chem. B, 1997, 101, 10506–10517.
34. E. S. Wiedner, J. Y. Yang, S. Chen, S. Raugei, W. G. Dougherty, W. S. Kassel, M. L. Helm, R. M. Bullock, M. R. DuBois and D. L. DuBois, Organometallics, 2012, 31, 144–156.
35. D. Devarajan and D. H. Ess, Inorg. Chem., 2012, 51, 6367–6375.
36. A. A. Franco, P. Schott, J. Jallut and B. Maschke, J. Electrochem. Soc., 2006, 153, A1053.
37. A. A. Franco, P. Schott, C. Jallut and B. Maschke, Fuel Cells, 2007, 7, 99.
38. K. Malek and A. A. Franco, J. Phys. Chem. B, 2011, 115, 8088–8101.
39. L. F. Lopes Oliveira, S. Laref, E. Mayousse, C. Jallut and A. A. Franco, Phys. Chem. Chem. Phys., 2012, 14, 10215–10224.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c3ra44896g

‡ Present address: 1) Laboratoire de Réactivité et de Chimie des Solides (LRCS), UMR CNRS 7314, Université de Picardie Jules Verne, 33 Rue Saint Leu, 80039 Amiens Cedex, France. 2) Réseau sur le Stockage Electrochimique de l’Energie (RS2E), FR CNRS 3459, France.

---