H2 activation by hydrogenase-inspired NiFe catalyst using frustrated Lewis pair: effect of buffer and halide ion in the heterolytic H-H bond cleavage.

Hydrogen is a clean fuel alternative to fossil fuels, and it is vital to develop catalysts for its efficient activation and production. We investigate the reaction mechanism of H2 activation in an aqueous solution by the recently developed NiFe complex (Ogo et al. Sci. Adv. 2020, 6, eaaz8181) using density functional theory (DFT) calculation. Our computational results showed that H2 is activated using frustrated Lewis pair. That is, H2 binds to the Fe site of the NiFe complex, acting as a Lewis acid, while the added buffer, acting as Lewis base, abstracts protons to form a hydride complex. Furthermore, the higher basicity in the proton abstraction reaction characterises reaction more exergonic and lowers the reaction barrier. In addition, in the proton abstraction by the water molecule, the reaction barrier was lowered when anion such as Cl- is in the vicinity of the water. Understanding the chemical species that contribute to the catalytic process in cooperation with the metal catalyst at the atomic level should help to maximise the function of the catalyst. This journal is © The Royal Society of Chemistry.

Introduction

Hydrogen is a clean energy fuel with zero carbon emissions from its combustion, catalytic systems that efficiently activate and generate hydrogen are in high demand. It is essential to build an environmentally friendly and sustainable catalytic system for the industrial application of catalysts. That is, catalysts are ideally created using earth-abundant elements that function in mild conditions with water as the solvent. Hydrogenase enzymes are known to activate/evolute hydrogen efficiently,[1-5] and their structure and function of hydrogenase have been investigated both experimentally and theoretically.[2,6] However, enzymes often have a stability issue. For example, [FeFe] hydrogenase is known to be unstable to the oxygen exposure.[7] Therefore, producing organometallic complexes that exceed the hydrogenase with respect to the robustness and efficiency is challenging.

Considering enzymatic reactions, the mechanism of hydrogen activation/evolution reaction by NiFe hydrogenase[8-10] is relatively less complicated than that of water splitting by photosystem II (PSII)[11-13] and ammonia synthesis from nitrogen by nitrogenase.[14,15] That is, because the active centre contains fewer metal elements than the active centre of PSII (Mn4CaO) and nitrogenase (FeMo-cofactor), the oxidation number fluctuates less, and the number of elementary reactions in the catalytic cycle of hydrogen activation and production is smaller. This could be related to the difficulty in designing/developing organometallic complexes. However, mimicking the function of H2 activation/generation and the catalytically active site using an organometallic catalyst is quite challenging in reality, and much effort has been put into developing a stable catalyst.[16-18]

The recently developed synthetic model for the catalytic centre of NiFe hydrogenase, [NiII(X)FeII(Cl)(CO)(L)](X = N,N′-diethyl-3,7-diazononane-1,9-dithiolate, L = 1,2,-bis(diphenylphosphino)ethane) (1, Fig. 1), reversely catalyses heterolytic H2 activation in aqueous solution[19] (Fig. 1). This catalyst is structurally much like hydrogenase in that a dithiolate ligand is used for bridging Ni and Fe, and the carbonyl ligand is coordinated to the Fe site. The H2 activation reaction efficiently occurs in the presence of phosphate buffer, and the hydride complex is obtained as a product (5, Fig. 1). X-ray crystallographic analysis for the hydride complex showed that the hydride ligand locates almost in the middle between Ni and Fe (5, Fig. 1) similar to the NiFe hydrogenase.[9] Although it was experimentally demonstrated that the NiFe complex has ability to catalyse H2 activation/evolution, the catalytic active site and the mechanism has not been clarified.

Fig. 1

H2 activation and evolution by NiFe complex.[19] A and B represent bases, and AH and BH represent conjugate acid of A and B, respectively.

Besides improving the catalyst itself, it is critical to investigate the effects of solvent, acid/base, and the other added chemicals on the catalytic reaction to develop and maximise the catalytic function. Alvarez-Hernandez et al.[20] reported that pKa of buffer acid affects the rate-determining step of the electrochemical hydrogenation reaction by Co-complexes in an aqueous solution. They demonstrated that the addition of buffer significantly increases the catalytic current, and the mechanism depends on pKa. By computationally testing several Brønsted acids with varied pKa values in the electrochemical CO2 reduction reaction with Re and Mn complexes, Riplinger et al.[21] demonstrated the effects of acids on thermodynamics and kinetics. They revealed that catalysis with weak acids required a more negative applied potential or higher acid concentration compared to catalysis with stronger acids.

Efficient transfer of protons is essential in hydrogen activation and evolution reactions. In this regard, it has been reported that incorporating pendant amines into the second coordination sphere of organic–metal complexes improves the efficiency of proton transfer.[22,23] While another strategy to enhance the efficiency of proton transfer should be to use water as a solvent that assists the proton relay. The importance of the water-assisted proton transfer is emphasised in the recent study for the electrochemical H2 production by cobalt-based complex.[24] The importance is not limited to the homogeneous catalyst. Our recent DFT study about alanine production from pyruvic acid at the TiO2 electrode also showed that the water molecule mediates the proton transfer from the TiO2 surface to the reactants[25] The construction of such an efficient catalyst system requires optimisation of many parameters, and it is critical to understand the influence of each parameter on the reaction free energy.

In this study, we first elucidate the mechanism of H2 activation by the NiFe complex. The reaction free energy and the reaction barrier are calculated using density functional theory (DFT) in conjunction with the implicit solvation model. The proton abstraction process is described by combining the implicit solvation model with explicitly added water molecules. In addition, we investigate the effects of buffer and halide ions in neighbouring water on the thermodynamics and kinetics in the step of proton abstraction from NiFe dihydrogen complex. We will also investigate the possibility that sulfur in the NiFe complex functions as a base. More specifically, we compare the following eight different bases; (1) [HPO4–(H2O)2]2−, (2) [H2PO4–(H2O)2]−, (3) [CH3COO–(H2O)]− (4) [H2O–(HPO4–H2O)]2−, (5) [H2O–(H2O)2], (6) [H2O–(Cl)]−, (7) sulfur in [NiFe]2+, and (8) sulfur in [NiFe–H2]2+. The species that hydrogen bonds with or interact electrostatically with the base are presented in parentheses. [NiFe]2+ in (7) and [NiFe–H2]2+ in (8) represent the dissociated complex of Cl− from complex 1 (Fig. 1), and NiFe dihydrogen complex, respectively. The proton transfer between NiFe complexes is considered in (7), while the proton transfer within the NiFe complex is considered in (8).

Method

All computations in this study were performed using the Gaussian 09 program.[26] All structures were fully optimised without any constraints using the BP86 functional[27,28] with Grimme's empirical dispersion corrections.[29] Our previous study has shown that the BP86 functional behaves well in this NiFe complex.[19] The diamagnetic properties of these ground states predicted by BP86 functional were consistent with 1H nuclear magnetic resonance spectra and the electron spin resonance. Further, the BP86 has also been used to study the mechanism of catalysis of NiFe hydrogenase, with reliable performance.[8,30] The Stuttgart/Dresden (SDD)[31] basis set and the associated effective core potential were employed for Ni and Fe, and def2-SVP basis sets[32] were used for the other atoms (BSI). The SMD implicit solvation model[33] was used to account for the solvation effects of water (ε = 78.4).

Vibrational frequency calculations were performed at the same level of theory to confirm the minima and TSs, and to obtain zero-point vibrational energy (ZPE) corrections. The thermal corrections were computed at 298.15 K and 1 atm pressure. Connectivity of the stationary points was confirmed by the ‘pseudo’ intrinsic reaction coordinate (IRC) approach,[34] where IRC calculations were performed for 20 steps from the TS (in both forward and backward directions), and subsequent structures were fully optimised to obtain the minima.

Potential energies of the optimised stationary points were calculated at the level of BP86-D3/BS2, (BS2 = SDD for Ni and Fe and def2-TZVP for the other toms). Solvation effects were considered using the SMD implicit solvation model with water as the solvent. Integrals were evaluated using the pruned grid consisting of 99 radial shells and 590 angular points per shell. The wave function stability was checked for all metal complexes.

Initial approximations of the TSs were obtained by the conventional AFIR methodology,[35] where the two-layer N-layered integrated molecular orbital and molecular mechanics (ONIOM) method was applied.[36] The ONIOM partitioning of the molecular system is shown in Fig. 2a. The BP86 functional was applied for the high level, using the SDD[37] basis set for Ni and Fe and the 3-21G basis sets[38] for the remaining atoms. The parameterisation method 6 (PM6)[39] was applied for the low level. In the AFIR simulation, an artificial force was applied between the atoms where bond formation is expected to occur (Fig. 2b). An artificial force parameter of 47.8 kcal mol−1 was used to explore the approximate reaction paths and TSs. The approximate TSs obtained by the AFIR method were finally optimised at the level of SMD/BP86-D3/BS1. The different use of basis set adapted in this study is that we employed previous studies[19,40] to avoid the high computational cost in the geometry optimisation and get the energy accuracy.

Fig. 2

(a) ONIOM partitioning; black and blue parts are treated at a high and low level, respectively. (b) Fragmentation in the AFIR study. The artificial force is applied between two red atoms.

Results and discussion

Analysis of complex 1 and complex 5

The experimentally determined oxidation states of complex 1 (Fig. 1) on Ni and Fe are both 2+, indicating that the orbital occupancies are d8 (Ni) and d6 (Fe), respectively. Therefore, to determine the ground state spin multiplicity of complex 1, we optimised the structures for the four spin multiplicities (S = 0–3), starting from the X-ray structure, and calculated the relative free energy (Table 1). Since the oxidation numbers of Ni and Fe for complex 5 (Fig. 1) are the same as those of complex 1, it is presumed that the complex 5 has the same orbital occupancy as complex 1.

Mulliken spin densities (ρ), enthalpies (ΔH)a, and free energies (ΔG)a for complexes 1 and 5 at four spin states

	S = 0	S = 1	S = 2	S = 3
Complex 1
ρ(Ni)	0.00	1.43	1.44	1.45
ρ(Fe)	0.00	0.00	1.94	3.45
ρ(S)	0.00	0.14	0.11	0.24
ρ(S)	0.00	0.14	0.16	0.27
ρ(Cl)	0.00	0.10	0.18	0.17
ΔH	0.0	16.0	39.4	59.2
ΔG	0.0	14.9	36.2	52.3
Complex 5
ρ(Ni)	0.01	1.27	1.23	1.35
ρ(Fe)	−0.01	0.25	2.11	3.29
ρ(S)	0.04	0.13	0.12	0.32
ρ(S)	−0.04	0.13	0.33	0.31
ρ(H)	0.00	0.02	0.00	0.16
ΔH	0.0	3.9	36.8	67.0
ΔG	0.0	3.3	34.8	61.9

Enthalpies and free energies are in kcal mol−1.

The calculated relative energies indicate that the ground state is singlet in both complexes 1 and 5. In complex 1, the singlet and triplet spin states are energetically well separated (∼15 kcal mol−1), but the energy difference is small in complex 5 (∼3 kcal mol−1). At higher spin multiplicities, S = 2 and 3, energy separation from the ground state is even greater. These results are consistent with previous calculations with acetone solvents.[19]

It is well known that it is difficult to correctly predict the energy difference between different spin states by DFT.[41] The spin multiplicity of the electron ground state is often mispredicted, especially when the energy split between spin states is small, such as in spin crossover complexes.[42] Therefore, it is important to investigate the density functional dependency for the system of interest. We investigated three local density functionals; M06-L,[43] BLYP[44,45]-D3, TPSS-D3,[3] and three hybrid density functionals; B3LYP[46]-D3 M06,[47] and TPSSh.[48] As a result, only BLYP-D3 correctly predicted that the ground state is a singlet (Table 2). TPSS-D3 shows that the singlet is the lowest energy as well as BLYP-D3. However, the energy difference between the singlet and triplet is rather small that it indicates that the two spin states are mixed, which is inconsistent with the experimental diamagnetic properties. Also, all hybrid functions incorrectly predict the ground state spin state.

Free energies (ΔG) of complex 5 for four spin states (S = 0, 1, 2, and 3) with different DFT functionals; BLYP-D3, M06-L, TPSS-D3, B3LYP-D3, M06, and TPSSh

	S = 0	S = 1	S = 2	S = 3
BLYP-D3	0.0	1.7	28.4	50.7
M06-L	0.0	−7.2	18.7	31.6
TPSS-D3	0.0	0.4	29.1	48.2
B3LYP-D3	0.0	−7.4	12.9	20.7
M06	0.0	−2.4	10.3	12.5
TPSSh	0.0	−4.6	16.8	26.6

Fig. 3 shows the optimised structures of complex 1 and 5 with selected structural parameters, and the more complete list of structural parameters are presented in Fig. S1.† Focusing on the distance between the metal and the ligand atom, the deviation from the X-ray structure for complex 1 is less than 0.1 Å for all selected bond lengths. On the other hand, in the optimized structure of the hydride complex 5, a deviation from the X-ray structure of 0.1 Å or more was observed between Ni and the hydride H. This large deviation in Ni–H bond length is due to the lack of chemical bonds, according to the Natural Bond Orbital (NBO) analysis. As shown in Fig. 4, the NBO with bonding orbitals between Fe and H is present, while the bonding orbital with electron density between Ni and H is absent. However, it is predicted that the hydride and Ni are electrostatically interacting with each other in complex 5 from the short NiFe distance in the hydride complex 5 (2.61 Å) compared to the complex 1 (3.23 Å). Such a reduction in metal–metal bond length due to electrostatic interactions between the Ni and Fe axial ligands is also seen in the O2 activation by the NiFe complex.[49] These interactions are important for intermediate stabilisation in the catalytic process.

Fig. 3

Key geometrical parameters of (a) complex 1 and (b) complex 5 in the optimised geometries.

Fig. 4

Natural bond orbital of NiFe hydride complex.

In our previously studied NiIr complex for H2 activation in which the ligand of Ni site is the same with the present NiFe complex,[50] X-ray structural analysis revealed that Cl− ligand coordinates to Ni. The coordination of Cl− or OH− stabilises the NiIr complex depending on the spin multiplicity, as subsequently confirmed by DFT calculation.[51] The coordination of OH− is unlikely because the H2 activation reaction by the present NiFe complex was conducted at an almost neutral pH. Also, the dissociated Cl− ligand rebinding in complex 1 is less likely due to the diluteness. Therefore, the binding of the solvent water molecule was considered.

Two positions, X and Y, were considered as the binding positions of the water molecules in complex 1 (Fig. 5). In the ground state, the binding to the positions X and Y are both impossible. Although the H2O binding state to the position X is located for S = 1–3, these states are thermodynamically unfavoured, as seen from the negative binding energy of H2O to the position X (Table 3). For the position Y, the stabilisation by the water binding is seen at S = 2 and S = 3, but these spin states are energetically much higher than the ground state (Table 1). These computational results eliminate the possibility of H2O binding for both X and Y positions.

Fig. 5

Two positions (X and Y) considered for H2O binding to complex 1.

Binding free energy (kcal mol−1) of a water molecule to the position X and Y in complex 1

	S = 0	S = 1	S = 2	S = 3
Position X
	—a	−9.9	−5.5	−5.0
Position Y
	—a	—a	1.0	2.7

H2O is unbound.

Mechanism of H2 activation by NiFe complex in aqueous solution

The H2 activation by NiFe complex was performed in the presence of phosphate buffer.[19] Because H2 activation was conducted at pH 7, eqn (2) is predominant among the three possible equilibria with respect to the protonation state of the phosphate buffer (eqn (1)–(3)). Therefore, the reaction mechanism was explored, assuming that HPO42− acts as a proton acceptor.

Fig. 6 shows the calculated free energy profile for H2 activation by NiFe complex. Starting from the complex 2 (Fig. 6), H2 binds to Fe in the η2 fashion via the transition state TS-23 (barrier height; 11.3 kcal mol−1). The buffer (HPO42−) then approaches the temporarily generated dihydrogen complex, 3, and abstracts a proton to form a hydride complex. This proton abstraction is an almost barrierless process.

Fig. 6

(a) Free energy profile of H2 activation by NiFe complex in the absence of buffer (HPO42−). ΔG and ΔH values (in parentheses) are given in kcal mol−1. (b) Optimized key geometries. The selected bond lengths are given in Å.

The free energy profile (Fig. 6) starts from water binding complex, 2, rather than complex 1. This description is more realistic than starting from complex 2 as it is less likely to recombine once Cl− is removed from the dilute complex 1. Another reason to exclude complex 1 is that the implicit solvation model does not quantitatively describe the solvation energy of the ionic species. According to our previous investigation for NiIr complex,[51] the hydration energy of Cl− was estimated as (−67.3 kcal mol−1) using SMD implicit solvation model, which is ∼6 kcal mol−1 lower than the experimental value (−73 ± 2 kcal mol−1).[52,53] As the free energy of complex 2 was calculated as 8.1 kcal mol−1 relative to complex 1 using the SMD model, the complex is estimated to be destabilised by Cl− dissociation by about 2 kcal mol−1, taking into account the error. Note that water molecules are unlikely to bind to the Fe site after Cl− dissociation, because water binding destabilises the complex by 8.0 kcal mol−1.

The binding site for H2 to the metal site is Fe, not Ni. This can be predicted from the molecular orbital diagram mainly including the 3d orbitals of Ni and Fe shown in Fig. 7. The electronic configuration is (Ni-3d)2(Ni-3d)2(Ni-3d)2(Ni-3d)2(Ni-3d)0 for Ni, and (Fe-3d)2(Fe-3d)2(Fe-3d)2 (Fe-3d)0 (Fe-3d)0 for Fe. The highest occupied orbit (HOMO) and the lowest empty orbit (LUMO) are on Ni and Fe, respectively, and there is no empty orbit at the vertical direction of Ni. In the optimised structure of dihydrogen complex (Fig. 6b), the H–H bond distance is 0.86 Å. This dihydrogen complex is classified as the most abundant Kubas-type complex[1] with an H–H distance of less than 1.0 Å and remaining H–H bonds.

Fig. 7

Molecular orbital diagrams (isovalue 0.02) of Ni and Fe for ligand-free NiFe complex at the vertical position, 2. The orbital energy is given in eV.

The H2 binding process has been shown to be an endergonic reaction (Fig. 6b), but it may be overestimated. Our previous study predicts that the main cause is an overestimation of translational entropy by treating H2 as an ideal gas.[51] Therefore, the reaction barrier of actual H2 addition is lower than calculated. Further, H–H bond cleavage is close to the barrierless process. Thus, it is predicted that H2 activation by NiFe complex is a very fast reaction at room temperature.

The binding energy of OH− to the Fe site was calculated to predict whether this catalyst works in alkaline conditions or not. The binding energy of OH− was estimated to be 13.3 kcal mol−1 against H2 binding energy (−8.6 kcal mol−1), indicating that the OH−coordination substantially stabilises the NiFe complex as opposed to the H2 binding. This binding energy comparison suggests that OH− inhibits H2 binding to Fe, which reduces efficiency under alkaline conditions. Such a large stabilisation of complex by binding of OH− to the metal centre was also found in the NiIr catalyst we theoretically investigated earlier,[51] and which is consistent with the experimental observation that the activity of the NiIr complex actually decreases under alkaline conditions.[50] The Ni–B state, which is a catalytically inactive state in NiFe hydrogenase, is recognized as the state having an OH− ligand bridging Ni and Fe.[54]

Comparison with H2 activation by NiFe–CH3CN and NiIr–Cl complexes

The revealed mechanism is similar to H2 activation by NiFe–CH3CN and NiIr–Cl complexes (Fig. 8a–c) previously investigated using DFT,[40,51] where ‘CH3CN’ in NiFe–CH3CN and ‘H2O’ in NiIr–H2O represent the coordination species in the active catalytic site Fe and Ir, respectively. In all three complexes, dihydrogen is activated using a frustrating Lewis pair. Common to all three catalyst systems, Lewis acids are metal complexes and Lewis bases are other additive reagents that do not react with each other (Fig. 8d–f). The dihydrogen weakly binds to a metal site other than Ni, then the H–H bond is heterolytically cleaved by the Lewis base.

Fig. 8

Two previously investigated dinuclear complexes catalyzing H2 activation (a) NiFe–CH3CN and (b) NiIr–Cl and current (c) complex 1. Square planar symmetry local to Ni commonly seen in all three complexes are circled by red dash lines. Optimized transition states of H–H bond cleavage by (d) NiFe–CH3CN,[40] (e) NiIr–Cl,[51] and (f) complex 1. Lewis acid and base are shadowed by green and red, respectively.

The mechanism similarity is expected from the structural analogy. All three complexes have a square planar symmetry local to Ni (Fig. 8a–c). Ni is chemically inert during the catalytic reaction, but it plays an important role in stabilising hydride complexes by electrostatically interacting with hydrides bound to Fe or Ir. In all three complexes, the H2 activation occurs with moderate strength of the base. Although strong base, CH3O−, was used in the experiment for H2 activation by the NiFe-CH3CN complex,[55] the subsequent computational studies have shown that H2 activation is also exergonic with the weak base CH3COO−.[56] Therefore, all three complexes tend to activate H2 rather than generate H2.

According to the ionization potential theorem, the negative of the HOMO energy corresponds to the ionization potential.[57] The calculated HOMO energy level of the ligand free complex (2, Fig. 7) is −4.87 eV which is higher than NiFe–CH3CN complex (−5.14 eV),[40] and NiIr–Cl complex (−5.27 eV),[51] suggesting that the present NiFe complex has more capability to donate an electron to the ligand. In NiFe–CH3CN and complex 1, phosphorus is coordinated to iron in both complexes, whereas the CO ligand is in complex 1 but not in NiFe–CH3CN. Therefore, the higher HOMO energy level in complex 1 than NiFe–CH3CN is probably due to this CO coordination. The CO ligand is also contained in NiFe and FeFe hydrogenase, and it maintains the low valence state of iron by the delocalization of pi electrons.[58]

Comparison with H2 activation by NiFe hydrogenase

Sigbhan et al.[59] recently used a cluster model involving active metal centre and several amino acid residues around it to theoretically study the mechanism of H2 activation by NiFe hydrogenase. One similar tendency for NiFe complexes and hydrogenases is the weak binding of H2 to the metal site. Another similarity is the use of frustrated Lewis pairs in the H–H bond cleavage. On the other hand, there are differences between hydrogenase and NiFe complexes. One is the active site, which is Fe in this catalyst but Ni in hydrogenase. The other distinction is the barrier to proton extraction. The NiFe complex abstracts protons without a reaction barrier. On the other hand, the barrier of H–H bond cleavage by NiFe hydrogenase is calculated to be ∼15 kcal mol−1 using the B3LYP functional, assuming that cysteine is the base.[59] The reaction barrier has been shown to depend on the DFT used and the percentage of Hartree–Fock exchange.[59] However, even considering them, it is expected that there is a substantial difference.

Effects of a base on thermodynamics and kinetics of H–H bond cleavage

Next, we investigate the effects of bases on the reaction energy of H2 activation and kinetic barrier in the proton abstraction step. The following eight reactions, (R1)–(R8), were examined.

The third term of each equation is the base for proton extraction, and the molecules in the vicinity of the species that directly extract the proton are shown in parentheses. The products and reactants of (R1) and (R4) are the same, but the species that directly extract protons in the transition state are different. (R7) is the reaction in which the sulfur of [NiFe]2+ (complex 2) abstracts protons of NiFe dihydrogen complex (3) to form hydride complex and (R8) is the reaction in which proton transfer occurs within the NiFe dihydrogen complex. The chemical structural formulas of the [NiFe (SH)]3+ in (R7) and the [NiFe(SH)–H]2+ in (R8) are shown in Fig. 9.

Fig. 9

Chemical structural formulas of (a) [NiFe(SH)]3+ in (R7) and (b) [NiFe(SH)–H]2+ in (R8).

The calculated free energy profiles for each H–H bond cleavage step in (R3)–(R6) are summarised in Fig. 10. The reaction energy, barrier height, and pKa of conjugate acid[60,61] of the employed base for (R1)–(R8) are summarised in Table 4. The barrier height was calculated as the free energy difference between the corresponding TS and complex 2, as the free energy level of intermediate complex 3 is theoretically overestimated as discussed in previous Subsection 3.2. A small pKa value for a conjugate acid indicates high acidity and low basicity of the base. According to a prior study, the DFT estimated pKa values deviate significantly from the experimental value. For example, the BP86 functional used in this study has been shown to give the mean unsigned error of 3.70 from experimentally observed pKa.[62] Moreover, the error is not systematic. Therefore, the discussion is based on the experimentally obtained pKa.

Fig. 10

Free energy profiles of proton abstraction by Lewis base from NiFe dihydrogen complex for four reactions; (a) for (R3), (b) for (R4), (c) for (R5), and (d) for (R6). ΔG and ΔH values (in parentheses) are in kcal mol−1. The selected atomic charges, q, are also shown.

Reaction free energy of H2 activation (kcal mol−1), barrier height of H–H bond cleavage (kcal mol−1), and experimental pKa value of the conjugated acid of the used Lewis base

	Reaction energy	Barrier heightc	Conjugate acid	pKa (exp.)d
(R1) [NiFe]2+ + H2 + [HPO4–(H2O)2]2− → [NiFe–H]+ + [H2PO4–(H2O)2]−a	−10.3	2.2	H2PO4−	7.2e
(R2) [NiFe]2+ + H2 + [H2PO4–(H2O)2]− → [NiFe–H]+ + [H3PO4–(H2O)2]−	11.2	—h	H3PO4	2.1e
(R3) [NiFe]2+ + H2 + [CH3COO–(H2O)]− → [NiFe–H]+ + [CH3COOH–(H2O)]	−2.4	9.4	CH3COOH	4.76f
(R4) [NiFe]2+ + H2 + [H2O–(HPO4–H2O)]2− → [NiFe–H]+ + [H2PO4–(H2O)2]−b	−10.3	8.8	—
(R5) [NiFe]2+ + H2 + [H2O–(H2O)2] → [NiFe–H]+ + [H3O–(H2O)2]+	4.8	16.8	H3O+	−1.74g
(R6) [NiFe]2+ + H2 + [H2O–(Cl)]− → [NiFe–H]+ + [H3O–(Cl)]	2.6	12.8	—
(R7) [NiFe]2+ + H2 + [NiFe]2+ → [NiFe–H]+ + [NiFe(SH)]3+	32.6	—h	—
(R8) [NiFe]2+ + H2 → [NiFe(SH)–H]2+	21.5	—h	—

HPO4− abstracts a proton from dihydride complex.

H2O abstracts a proton from dihydride complex.

Free energy difference between the corresponding TS and complex 2.

pKa value of the conjugate acid.

The value taken from ref. 64.

The value taken from ref. 63.

The value taken from ref. 66.

Transition state is not determined.

The reaction (R1) is already discussed in Section 3.2. Then, we next investigated whether H2PO4−, which coexists with HPO42− in aqueous solution (eqn (2)), functions as a base (R2). The calculated reaction free energy for (R2) is endergonic by 11.2 kcal mol−1 (Table 4) and the reaction is less likely to occur.

The case of CH3OO− acting as a proton acceptor (R3) is shown in Fig. 10a. The conjugated acid of CH3COO−, CH3COOH, is smaller pKa (pKa = 4.8)[63] than H2PO4− (pKa = 7.2),[64] indicating that HPO4− shows the larger basicity than CH3COO−. Compared to the case where HPO4− acts as a base (−10.3 kcal mol−1), the reaction is less exergonic (−2.4 kcal mol−1). Furthermore, the barrier height is higher than when HPO42− is used as a base. The same tendency is seen in previously studied NiFe complex. The proton abstraction step is more exergonic when CH3O− is used (∼–30 kcal mol−1)[40,56] as the base than when CH3COO− (∼–10 kcal mol−1)[56] is used. Similar trends have been reported[65] in the activation of hydrogen molecules by a metal-free catalyst, where 75 Lewis pairs have been tested using DFT calculations.[65]

Next, we examined the reaction (R4) in which H2O abstracts protons instead of HPO42− (Fig. 10b). In this model reaction, the transition state was calculated assuming that HPO42− is hydrogen-bonded to the water molecule (Fig. 10b). Note that the final product is H2PO4− + H2O rather than HPO42− + H3O+, as the former is thermodynamically advantageous. The calculated barrier height is 8.8 kcal mol−1, which is higher than the direct proton abstraction by HPO42− (2.2 kcal mol−1). This can be inferred from the smaller pKa (−1.74)[66] of H3O+ than that of H2PO4−. As a result, the abstraction by H2O is kinetically unfavoured compared to proton abstraction by HPO42−. However, the water molecule has more opportunities to be placed near the metal complex than HPO4− due to the concentration, and it would be possible that the water molecule abstracts proton.

Only water can act as a base in the absence of buffer (R5). Therefore, we added two explicit water molecules in addition to the water molecules that abstract the proton in the model reaction (Fig. 10c). The barrier is higher (16.8 kcal mol−1) than when HPO42− (2.2 kcal mol−1) or CH3COO− (9.4 kcal mol−1) extracts proton. This reaction barrier can be overcome, but the reaction is less likely to occur as the it is endergonic.

In the model reaction of (R6), the Cl− halide ion locates the vicinity of the water molecule, which abstracts the proton (Fig. 10d). The proximity of Cl− to the water substantially lowers the reaction barrier (∼4 kcal mol−1), demonstrating the kinetic advantage. In the transition state, the Mulliken charge of oxygen in the water that abstracts the protons is confirmed to be more negative when Cl− (Fig. 10d) is near the water. From this, it is expected that the proton affinity of water will increase when negatively charged ions are present near water molecules. However, this process is an endergonic reaction, and it is expected that no apparent reaction will occur.

In hydrogenase, it has been proposed that cysteine sulfur functions as a proton receptor.[2,67] Since sulfur that crosslinks Ni and Fe is involved in the NiFe complex, we investigated whether it functions as a base ((R7) and (R8)). As a result, the calculated reaction free energy was substantially an endergonic reaction (Table 4) for both when proton transfer occurred between the NiFe complexes (R7) and when proton transfer occurred within the NiFe complex (R8). Therefore, the sulfur in the NiFe complex does not function as a proton acceptor.

Overall, both the reaction energy of heterolytic H–H bond cleavage and the barrier height simply correlate with the strength of the proton abstracting base. Therefore, the use of stronger bases is beneficial for H2 activation. It was also found that the reaction barrier is also affected by anionic species near the water molecules that extract protons. On the other hand, it has been shown that bases such as OH− may interfere with H2 activation by binding strongly to the active site, suggesting that base selection is not straightforward. More ideal is the H2 activation catalytic reaction, in which solvent water molecule acts as a base. Thermodynamic stability of the hydride complex is required to overcome the destabilisation of H3O+ formation to achieve an ideal H2 activation catalytic reaction. Another approach would be to reduce or oxidise the hydride complex irreversibly.

H2 evolution reaction

Fig. S2† shows the free energy profile of hydrogen production using acetic acid as the proton donor. Starting with the hydride complex, the proton transfer from acetic acid to Fe occurs first, then the dihydrogen complex is formed and H2 is released from Fe. Finally, the conformational change of the ligand (2 → 2′) stabilizes the complex by approximately 2 kcal mol−1. The complexes 2 and 2′ were obtained by structurally optimizing the complex in which the axial ligand was removed from the complex 1 and the complex 5, respectively.

DFT calculations show that the reactant (hydride complex) and the product (complex 2′) is energetically close and the reaction reversibly occurs. However, if the reaction begins with a hydride complex, the maximum amount of dihydrogen produced is the same concentration as the hydride complex utilized. Therefore, in reality, H2 activation is less likely to occur due to the low H2 concentration.

In contrast to the hydrogen activation reaction, stronger acids are more stabilized by proton abstraction and are thermodynamically preferred to facilitate the hydrogen production reaction. In fact, Ahmed et al.[68] showed that in the hydrogenation reaction with the developed NiFe complex, catalytic activity appears by adjusting to a lower pH.

Conclusions

DFT calculation was employed to elucidate the reaction mechanism of H2 activation by NiFe complex. The heterolytic H2 activation involves two steps: (1) formation of a dihydrogen complex by the coordination of H2 to the Fe site and (2) formation of a hydride complex by the abstraction of proton with a base.

The similarity with NiFe hydrogenase is that the frustrated Lewis pair is used for H–H bond cleavage. The Ni moiety certainly functions as a Lewis acid in hydrogenase, although arginine,[10,69] glutamine,[70,71] and cysteine,[2,67] are listed as potential candidates for Lewis bases, but have yet to be validated. Another similarity is that dihydrogen binds very weakly to the metal centre in the formation of dihydrogen complexes. As a result, avoiding extreme stabilization of intermediates is generally essential in the catalytic process. On the other hand, the difference is that Fe, rather than Ni, is the active site throughout the catalytic cycle, unlike NiFe hydrogenase. For the present NiFe complex, it was analysed that Ni is inert through the catalytic cycle because of its square planar symmetry and the absence of the empty orbital for the axial direction. In recent years, NiFe hydrogenase-inspired models for Ni-based proton reduction reactions have been proposed.[72] In this way, creating an artificial model that is closer to the enzyme enables a close comparison with the function of hydrogenase. It is necessary to change the ligand having the square symmetry of Ni to improve the present NiFe complex to a Ni-based catalyst.

The higher basicity of the Lewis base lowers the reaction barrier and makes the reaction more exergonic in the proton abstraction reaction from a dihydrogen complex. On the other hand, it is thermodynamically advantageous to use a stronger acid in the hydrogen generation reaction. We also found that when water abstracts protons, the reaction barrier becomes lower with the present of buffer bases or Cl− near the water. The contribution of buffered bases to catalytic reactions has also been reported in the formation of O–O bonds in water splitting.[73] In addition, our previous study has shown that the added ionic species participates in the catalytic process and lowers the reaction barrier of CO bond cleavage of carbon dioxide.[74] Thus, in addition catalysts, other chemical species used in catalytic reactions are often dominant in the catalytic process, and a systematic understanding at the atomic level is essential.

The experimental study for the present NiFe complex has shown that hydride complex generated by H2 activation undergoes isomerization reactions, and that these isomers have different reactivity in electron transfer, hydride transfer, and hydrogen generation reactions.[19] In the future, we will establish the mechanism of the isomerization reaction and elucidate whether the reactivity for the hydride transfer and H2 evolution is controlled by kinetics or thermodynamics.

Conflicts of interest

There are no conflicts of interest to declare.