H2 binding to the active site of [NiFe] hydrogenase studied by multiconfigurational and coupled-cluster methods

[NiFe] hydrogenases catalyse the reversible conversion of molecular hydrogen to protons and electrons. This seemingly simple reaction has attracted much attention because of the prospective use of H2 as a clean fuel. In this paper, we have studied how H2 binds to the active site of this enzyme. Combined quantum mechanical and molecular mechanics (QM/MM) optimisation was performed to obtain the geometries, using both the TPSS and B3LYP density-functional theory (DFT) methods and considering both the singlet and triplet states of the Ni(II) ion. To get more accurate energies and obtain a detailed account of the surroundings, we performed calculations with 819 atoms in the QM region. Moreover, coupled-cluster calculations with singles, doubles, and perturbatively treated triples (CCSD(T)) and cumulant-approximated second-order perturbation theory based on the density-matrix renormalisation group (DMRG-CASPT2) were carried out using three models to decide which DFT methods give the most accurate structures and energies. Our calculations show that H2 binding to Ni in the singlet state is the most favourable by at least 47 kJ mol . In addition, the TPSS functional gives more accurate energies than B3LYP for this system.


Introduction
Hydrogenases are metalloenzymes that catalyse the reversible cleavage of H 2 to protons and electrons. They are widespread in nature and found in bacteria, archaea, and some eukaryotes. 1 They can be classified into [NiFe], [FeFe], and [Fe] hydrogenases based on the metal-ion content in the active sites. The binuclear active site of [NiFe] hydrogenases is shown in Fig. 1. The iron ion is coordinated by one carbon monoxide and two cyanide molecules. Moreover, two thiolates from Cys-84 and 549 (the residues are numbered according to the enzyme from Desulfovibrio vulgaris Miyazaki F 2 ) bridge the two metals. Two additional cysteine ligands (Cys-81 and 546) coordinate terminally to the nickel ion. This gives a penta-coordinate iron ion and a tetracoordinate nickel ion, which are ready to bind an additional ligand during the reaction mechanism (H 2 or H À ) or as an inhibitor (e.g. OH À or CO). Typically, the [NiFe] hydrogenases also contain three iron-sulfur clusters and an octahedral Mg 2+ site.
The [NiFe] hydrogenases have been extensively studied by spectroscopic, 1,3,4 electrochemical, 5 biomimetic, 6,7 and computational methods. [8][9][10][11][12][13][14] Based on these studies, a tentative reaction mechanism has been suggested, shown in Fig. 2. The Fe ion is supposed to remain in the low-spin +II state throughout the reaction. In the Ni-SI a state, the Ni ion is in the +II oxidation state, without any extra ligand. Then, one electron and one proton are added to form the Ni-L state, in which the proton probably is bound to the S atom of the terminal Cys-546 residue. Next, the proton takes up two electrons from the Ni ion to form a hydride ion and moves to a bridging position between the two metals, thereby forming the Ni-C state. After that, another H + /e À pair is added, resulting in the Ni-R state, in which the proton again goes to Cys-546. Next, a H-H bond is formed between the hydride and the proton, giving a H 2 complex. Finally, the hydrogen molecule is released, regenerating the Ni-SI a state, which is ready to start a new reaction cycle.
In this mechanism, one of the key questions is how H 2 binds to the active site. Previous studies have suggested two different binding modes (Fig. 2), viz. binding to Ni or binding to Fe. Moreover, the ground state of the Ni ion is not clear when H 2 binds to the active site. Experimentally, the evidence is varying. Carbon monoxide, which is a competitive inhibitor of H 2 , binds to Ni. 15 Likewise, a likely binding path for H 2 has been identified by xenon binding experiments and it ends at the Ni ion. 16,17 On the other hand, from an organometallic perspective, Fe is the expected binding site of H 2 . 18 This would also explain the unusual choice of Fe ligands (CO and CN), which enforces a low-spin state on Fe and enhances H 2 activation. 19 Likewise, previous theoretical studies have given varying and often conflicting results, as is described in Table 1. In the first QM-cluster studies with minimal models, only H 2 binding to Fe in the singlet state was reported. 20,21 QM-cluster calculations with the density-functional theory (DFT) B3LYP method and somewhat larger models gave the same results and no local minimum for H 2 bound to Ni was found. 10 However, in the transition state for H-H cleavage, H 2 moved to mainly bind to Ni, but with one atom bridging to Fe (Ni-H distance of 1.68 and 1.77 Å and a Fe-H distance of 1.80 Å). Wu and Hall made a systematic study of the H 2 binding to QM-cluster models of varying sizes. They also pointed out that Ni II can attain either a low-spin singlet state or a high-spin triplet state. 22 For their smallest model (42 atoms), they found only Fe binding. However, with a larger model (103 atoms), they could find both binding modes in the singlet state. It turned out that H 2 preferred to bind to Fe in the triplet state and to Ni in the singlet state. However, the triplet state was always more stable by 5-52 kJ mol À1 . If Ni was oxidised to low-spin Ni(III), both binding modes were found, but binding to Fe was preferred by 1-21 kJ mol À1 .
Hillier and coworkers performed both QM-cluster and QM/MM calculations with a minimal QM system in both spin states and employing two DFT functionals BP86 and B3LYP. 23 The results were somewhat varying, but the QM-cluster calculations indicated that Fe-binding mode was most favourable, whereas Ni binding was preferred in most QM/MM calculations (except with B3LYP in the triplet state). Moreover, in the Fe-binding mode in the QM/MM structures, one of the H atoms bridged to Ni (1.82-2.11 Å distance). In all B3LYP calculations, the triplet state was most stable, whereas the opposite was true for the BP86 functional.
Recently, H 2 binding to the Ni-SI state was studied by Bruschi et al. with the BP86 functional using models of three   Two years ago, we studied the geometries and singlet-triplet energy splitting for the Ni-SI a state 25 using coupled-cluster calculations with singles, doubles and (perturbatively treated) triples (CCSD(T)), 26,27 multiconfigurational second-order perturbation theory with a complete active-space wavefunction (CASPT2) 28 as well as restricted active-space theory (RASPT2). 29 We also benchmarked various DFT methods to describe the geometries and singlet-triplet differences. TPSS was found to be a proper method for this system. In this paper, CCSD(T) and density matrix renormalisation group (DMRG) CASPT2 calculations [30][31][32] were performed to study the H 2 binding in the active site of [NiFe] hydrogenase based on structures obtained by QM/MM optimisation. We also used big-QM 33 calculations to estimate the influence of the surroundings.

Models
The calculations were started from our previous QM/MM structure of the Ni-SI a state of [NiFe] hydrogenase. 14 It was based on the 1.4 Å crystal structure from D. vulgaris Miyazaki F (PDB code 1H2R). 2 The structure of the Ni-SI a state is shown in Fig. 1. H 2 was then added to the active site, giving a total of 56 989 atoms. The protonation states of all residues were the same as in the previous structure of the deprotonated Ni-SI a state. 14 In this paper, two binding modes of H 2 were examined, viz. binding to Ni or to Fe (Fig. 2).
For the DFT and DMRG-CASPT2 calculations, we used models of three different sizes, shown in Fig. 3, whereas for the CCSD(T) calculations, only the minimal model 1 could be afforded. Model 3 (47 atoms) contained the Ni and Fe ions with their first sphere ligands (CO, two CN À , and four Cys residues, modelled by CH 3 S À ), an acetic-acid model of Glu-34, an imidazole model of His-88, and the H 2 molecule. It was taken from structures optimised with QM/MM. Model 2 (30 atoms) was obtained from these structures by deleting the acetic-acid and imidazole groups, whereas model 1 (18 atoms) was obtained by replacing the methyl groups by a proton, employing a S-H bond length of 1.34 Å. Thus, only single-point calculations in vacuum were performed.

QM calculations
DFT calculations were performed with the Turbomole 7.0 software. 34 We employed two DFT methods, TPSS 35 and B3LYP [36][37][38] with the def2-QZVPD 39,40 basis set. They represent typical and widely used examples of one pure and one hybrid DFT functionals. The calculations were sped up by expanding the Coulomb interactions in an auxiliary basis set, the resolution-of-identity (RI) approximation. 41,42 Empirical dispersion corrections were included with the DFT-D3 approach, standard zero-damping, and optimised parameters for each DFT method, 43,44 as implemented in Turbomole.

View Article Online
The triplet state was obtained with UHF formalism. Spin contamination was small (hS 2 i = 2.007-2.014). DMRG-CASPT2 calculations [30][31][32] were performed with BLOCK code 45 54 The number of renormalized states (m) was 1000 in all calculations. We used the default sweep schedule implemented in BLOCK, i.e. the initial DMRG iterations were done with a small m value to approximately converge the DMRG wavefunction, whereas the later iterations were carried out with the full m. A small amount of perturbative noise (e = 10 À4 in BLOCK) was added to the wavefunction in the initial iterations to prevent it from being trapped in a local minimum. The orbital ordering was automatized by minimising the quantum entanglement using a genetic algorithm. 55 In all DMRG-CASPT2 calculations, the four-particle reduced density matrices were obtained by a cumulant reconstruction based on the two-and three-particle reduced density matrices, neglecting the four-cumulant (the DMRG-cu(4)-CASPT2 approximation; this approach will be called DMRG-CASPT2 throughout this article). 31 DMRG-cu(4)-CASPT2 calculations with m = 1000 have been shown to work well for several transitionmetal models with medium-sized active spaces, including [NiFe] hydrogenase models. 56 All DMRG-CASPT2 calculations were performed with the standard ionization-potential electron-affinity (IPEA) Hamiltonian shift of 0.25 a.u. 57 An imaginary level shift of 0.1 a.u. 58 was used to avoid weak intruder states and to improve convergence of the perturbational treatment. The core electrons of the C, N, O, and S atoms were kept frozen in the perturbational treatment, whereas all valence electrons were correlated. To investigate the impact of the Fe and Ni semi-core (3s,3p) correlation, we performed calculations both with and without the 3s and 3p electrons correlated. Scalar relativistic effects were included using a standard second-order Douglas-Kroll-Hess Hamiltonian.
Spin-restricted CCSD(T) calculations 26,27 were carried out with the MOLPRO 2012.1 software. 59 These calculations employed the same basis set as the DMRG calculations. We also performed calculations both with and without the 3s and 3p electrons correlated for the Ni and Fe ion. The CCSD(T) calculations are too demanding to be used with model 2 and 3. Therefore, the CCSD(T) energies presented in this paper for models 2 and 3 were obtained by adding the difference between the CCSD(T) and DMRG-CASPT2 result for model 1 to the DMRG-CASPT2 energies (without the 3s, 3p semi-core correlation) for these models: We will see below that the correction (last two terms in eqn (1)) is rather small, 0-13 kJ mol À1 .
Thermal corrections to the Gibbs free energy (DDG therm ; including the zero-point vibrational energy) were calculated at 298.15 K and 1 atm pressure using an ideal-gas rigid-rotor harmonic-oscillator approach 60 from vibrational frequencies calculated for model 3 in vacuum at the TPSS/def2-SV(P) level of theory. To obtain more stable results, low-lying vibrational modes (below 100 cm À1 ) were treated by a free-rotor approximation, using the interpolation model suggested by Grimme and implemented in the thermo program. 61 QM/MM calculations QM/MM calculations were performed with the ComQum software. 62,63 In this approach, the protein and solvent are split into two subsystems: system 1 (the QM region) was relaxed by QM methods and it consisted of model 3 in Fig. 3. System 2 consisted of all the remaining parts of the protein and the solvent (577 crystal-water and 14 311 bulk-water molecules). It was kept fixed at the equilibrated crystal structure. Our previous study of [NiFe] hydrogenase showed that the dynamic effect of the surrounding protein (obtained by calculating QM/MM free energies) is restricted (up to 21 kJ mol À1 ). 25 In the QM calculations, system 1 was represented by a wavefunction, whereas all the other atoms were represented by an array of partial point charge, one for each atom, taken from MM libraries. Thereby, the polarization of the QM system by the surrounding protein is included in a self-consistent manner (electronic embedding). When there is a bond between system 1 and 2 (a junction), the hydrogen link-atom approach was employed: the QM system was capped with hydrogen atoms (hydrogen link atoms, HL), the positions of which are linearly related to the corresponding carbon atoms (carbon link atom, CL) in the full system. 62,64 All atoms were included in the pointcharge model, except the CL atoms. 65 The total QM/MM energy in ComQum was calculated as 62,63 where E HL QM1+ptch2 is the QM energy of the QM system truncated by HL atoms and embedded in the set of point charges modelling systems 2 (but excluding the self-energy of the point charges). E HL MM1;q 1 ¼0 is the MM energy of the QM system, still truncated by HL atoms, but without any electrostatic interactions. Finally, E CL MM12;q 1 ¼0 is the classical energy of all atoms in the system with CL atoms and with the charges of the QM system set to zero (to avoid double counting of the electrostatic interactions). By this approach, which is similar to the one used in the ONIOM method, 66 error caused by truncation of the QM system should cancel.
The geometry optimisation was continued until the energy change between two iterations was less than 2.6 J mol À1 (10 À6 a.u.) and the maximum norm of the Cartesian gradients was below 10 À3 a.u. The QM calculations were performed with Turbomole 7.0 software. 34 Geometry optimisations were performed with the TPSS 35 and B3LYP 36-38 functionals in combination with def2-SV(P) 67 basis set, including empirical dispersion corrections with the DFT-D3 approach. 68 The MM calculations were performed with Amber software, 69 using the Amber ff14SB force field. 70

Big-QM calculations
Previous studies of [NiFe] hydrogenase have shown that both QM-cluster and QM/MM energies strongly depend on the size of QM system. 13,65,71 To avoid this problem, we have developed the big-QM approach to obtain converged results: 33 the QM system was extended with all residues with at least one atom within 4.5 Å of model 2 in Fig. 3 and junctions were moved at least two residues away. In addition, all charged groups buried inside the protein were included, but the three iron-sulfur clusters were omitted to avoid convergence problems, which according to our previous calculations can be done without compromising the energies. 33,72 This gave a QM system of 819 atoms, shown in Fig. 4. 14 All the big-QM calculations were performed on coordinates from the QM/MM calculations and with a point-charge model of surroundings, because this gave the fastest calculations in our previous tests. 33 They also employed the multipole-accelerated resolution-of-identity J approach (marij keyword). Two different functionals, TPSS and B3LYP, with def2-SV(P) basis set, were used (the same as in the QM/MM optimisations).
To this big-QM energy, we added the DFT-D3 dispersion correction, calculated for the same big QM system with Becke-Johnson damping, 44 third-order terms, and default parameters for the TPSS or B3LYP functional using dftd3 program. 43,44,68 We also included a standard QM/MM correction for this large QM system (eqn (2), but with the big-QM region).
Reported final energies are the big-QM energies, including dispersion and the MM correction (E TPSS/SV(P) bigQM/MM ). This energy was extrapolated to the CCSD(T) level using QM calculations on model 3 in gas phase and thermal corrections were added:

Results and discussion
In this paper, we study the binding of H 2 to the active site of [NiFe] hydrogenase, addressing three questions: Does H 2 prefer to bind to Fe or Ni? Does the binding take place in the singlet or the triplet state? Which DFT functional gives the more reliable energies? This is done by employing two accurate QM methods, CCSD(T) and DMRG-CASPT2. We start by discussing geometries for the structures obtained with two DFT functionals, TPSS and B3LYP. Then, we analyse the orbitals from DMRG calculations. Finally, the relative energies of different binding modes and spin states will be compared.

Geometries
Structures were first optimised by QM/MM, using the 47-atom QM system of model 3 in Fig. 3. The QM/MM optimisations were performed with the TPSS and B3LYP functionals using the def2-SV(P) basis set. For the singlet (low-spin) state, two binding modes of H 2 were located (shown in Fig. 5     Several precautions have been taken to ensure that the presented structures are accurate and reliable. First, we reoptimised the B3LYP structures starting from the corresponding TPSS structures and vice versa to ensure that the reported differences are significant. Second, the structures were reoptimised by QM/MM calculations, in which all residues with at least one atom within 6 Å of any atom in the QM region were relaxed by MM. This did not change the structures significantly, and the same binding modes were still located. Third, the observation that only the Fe-H 2 binding mode is found for the triplet state is in accordance with all previous studies, except the QM-cluster calculations with a 290-atom model, for which also the 3 Ni-H 2 state was reported. 13 Therefore, we performed QM/MM optimisations, using a 298-atom QM region (we use [NiFe] hydrogenase from a different species, so a few residues are different). However, we could still not locate the 3 Ni-H 2 state, even if the Ni-H A and Ni-H B bonds were initially restrained to the distances in their study. On the other hand, we could find the 1 Fe-H 2 state with the same QM/MM method and the 298-atom QM region, although this state was not found in the previous study.

Orbital analysis in the DMRG calculations
In the DMRG calculations, we used the three models shown in Fig. 3. In all three models, the same active space was employed. For the singlet state, the active space contained the two H 2 1s orbitals, five Ni 3d orbitals, five Fe 3d orbitals, and a second set of correlating d orbitals for each occupied metal orbital

View Article Online
(which is often referred to as 3d 0 ) to account for the so-called double-shell effect (four on Ni and three on Fe). 74 In addition, three bonding orbitals between Ni, Fe, and their ligands were included. Thus, the active space consisted of 22 electrons in 22 orbitals. The active orbitals from the DMRG calculation are shown in Fig. 6 for the 1 Ni-H 2 state, obtained from the QM/MM optimisation with the TPSS functional. Three doubly occupied orbitals (orbitals 1, 3, and 8) show interactions between H 2 and the Ni ion. The first two orbitals involve the H-H s bonding orbital interacting with a Ni 3d orbital and ligand bonding orbital, whereas the third involves the H-H s* antibonding orbital interacting with another Ni 3d orbital. There is also a weak interaction between the H-H s* antibonding orbital and the Fe 3d orbital in orbital 8. Similar orbitals were obtained for the same state with the B3LYP structures (Fig. S1, ESI †).
Next, we studied the 1 Fe-H 2 state obtained from the QM/MM optimisations with both the TPSS and B3LYP functionals. The orbitals for the former are collected in Fig. 7, whereas the orbitals of the latter are shown in Fig. S2 (ESI †). It can be seen from Fig. 7 that the H-H s bonding orbital interacts with Ni in orbitals 7, with Fe in orbital 3, and with both metals in orbital 1.
For the triplet 3 Fe-H 2 state, four 3d orbitals can be involved in metal-ligand bonding (two singly occupied Ni 3d orbitals and two empty Fe 3d orbitals), whereas in the singlet state there are only three such orbitals (one empty Ni 3d orbital and two empty Fe 3d orbitals). Therefore, an extra ligand orbital was included in the active space. Moreover, as all five Ni 3d are occupied in the triplet state (instead of four in the singlet states), the active space was also extended with an extra 3d 0 orbital. This gives an active space of 24 electrons in 24 orbitals (Fig. 8). The H-H s bonding orbital interacts only weakly with Ni (orbital 1). On the other hand, there is a strong interaction between H 2 and Fe in orbital 6 and a weaker interaction also in orbital 14.

Energies
Next, we focus on the energies. We have studied the three states ( 1 Ni-H 2 , 1 Fe-H 2 and 3 Fe-H 2 ) using QM/MM geometries obtained with either the TPSS or the B3LYP method with the def2-SV(P) basis set. All energies discussed come from singlepoint energies calculated on these six QM/MM structures. Energies were calculated with four different methods: TPSS-D3 and B3LYP-D3 with the much larger def2-QZVPD basis set, as well as CCSD(T) and DMRG-CASPT2 with the even larger ANO-RCC basis set. The calculations were performed on the three models in Fig. 2 in vacuum (CCSD(T) calculations could be afforded only for model 1). The relative energies of the six states are shown in Table 3, always using the 1 Ni-H 2 state as the reference. In the CCSD(T) calculations, the 3s,3p correlation effect of the metal atoms was included. In a recent study, we have demonstrated that CASPT2 often fails to correctly describe the contribution of 3s,3p correlation energetics. 75 Therefore, in the DMRG-CASPT2 results shown in Table 3, the 3s,3p correlation was not obtained with perturbation theory, but it was taken from CCSD(T) calculations with and without 3s,3p correlation. A table comparing the CASPT2 and CCSD(T) relative energies obtained with or without metal 3s,3p correlation is provided in the ESI † (Table S1). For the present relative energies, the effect of the 3s,3p correlation was quite small, up to 13 kJ mol À1 , and the difference in 3s,3p correlation between CCSD(T) and DMRG-CASPT2 was 4-7 kJ mol À1 .
We first discuss results obtained with model 1, i.e. the smallest model, in which the cysteine residues are modelled by SH À . For this model, we could obtain energies with all four methods, TPSS, B3LYP, CCSD(T), and DMRG-CASPT2. The occupation numbers in Fig. 6-8 and Fig. S1-S3 (ESI †) show that all three studied states are essentially single configurational (they differ by less 0.09 from 0, 1, or 2). This is also supported by the T 1 , D 1 , and %TAE diagnostics 76 for the CCSD(T) calculations, shown in Table S2 (ESI †). They are similar to our previous CCSD(T) calculations on bimetallic [NiFe] hydrogenase models. 25 Therefore, we will use the CCSD(T) results as the reference.
From Table 3, it can be seen that TPSS predicts 1 Ni-H 2 as the most stable binding mode, but only 3 kJ mol À1 more stable than the 3 Fe-H 2 state and 26 kJ mol À1 more stable than the 1 Fe-H 2 state. The TPSS results are in satisfactory agreement with CCSD(T) in predicting the most stable binding mode. DMRG-CASPT2 overestimates the stability of the 3 Fe-H 2 state by 13 kJ mol À1 compared to CCSD(T). On the other hand, B3LYP gives quite different results, strongly favouring the triplet 3 Fe-H 2 state (by 55 kJ mol À1 for the B3LYP structure). For the B3LYP structures, CCSD(T) and DMRG-CASPT2 partially fix this problem by significantly destabilising the triplet state compared to B3LYP, reducing the energy difference between 1 Ni-H 2 and 3 Fe-H 2 by B40 kJ mol À1 . However, this destabilisation is not large enough to bring the 1 Ni-H 2 state back as the ground state. On the other hand, the CCSD(T) and DMRG-CASPT2 energy results agree within 2 kJ mol À1 for the B3LYP structures. Thus, the two ab initio methods give quite similar results, but it is not clear whether 1 Ni-H 2 or 3 Fe-H 2 is the most stable structure. Fortunately, this ambiguity disappears for the larger and more realistic models.
For models 2 and 3, no CCSD(T) calculations could be afforded. Instead, the CCSD(T) results in Table 3 were simply obtained by correcting the DMRG-CASPT2 results with the difference of the CCSD(T) and DMRG-CASPT2 results for model 1 (eqn (1)). From Table 3, it can be seen that the 3 Fe-H 2 state is strongly destabilised with the larger models. It is now the least stable state, 26-46 kJ mol À1 less stable than the 1 Ni-H 2 state, according to the combined CCSD(T) and DMRG-CASPT2 energies. B3LYP still favours the triplet state, but these energies are wrong by 58-69 kJ mol À1 , according to the ab initio data.
The DMRG-CASPT2 energies change by up to 22 kJ mol À1 when model 2 is increased to model 3 (i.e. adding models of two second-sphere residues, hydrogen-bonding to the active site), showing that the energies strongly depend on the surrounding protein. To more accurately model the surroundings, we have used the big-QM approach. Our big-QM model included 819 atoms (Fig. 4), which is much more than in previous studies of H 2 binding to [NiFe] hydrogenase (19-290 atoms, cf. Table 1). It can be seen from Table 4 that the big-QM energies (E TPSS/SV(P) bigQM/MM ) favour the 1 Ni-H 2 state by 50 and 53 kJ mol À1 compared to the 1 Fe-H 2 and 3 Fe-H 2 states. These large energy differences come mainly from the QM energy contribution: the MM corrections (E CL;big QM MM12;q 1 ¼0 À E HL;big QM MM1;q 1 ¼0 ) are small, 2-3 kJ mol À1 (owing to the large big-QM region), whereas the dispersion contributions are 9-11 kJ mol À1 . This energy (obtained at the TPSS/def2-SV(P) level) should be corrected by the CCSD(T) and DMRG-CASPT2 results for models 1 and 3. These corrections are slightly larger than in Table 3 because they start from TPSS calculations with the def2-SV(P) basis set. Finally, thermal corrections to Gibbs free energy were added (DDG therm ), obtained from frequency calculations on the isolated model 3. They favour the 3 Fe-H 2 state by 7 kJ mol À1 . With all corrections added (G tot in Table 4), 1 Ni-H 2 is 47 and 78 kJ mol À1 more stable than the other two states.

Conclusions
In this paper, we have studied the binding of H 2 to the active site of [NiFe] hydrogenase with coupled-cluster CCSD(T) and DMRG-CASPT2 methods. The geometries were optimised by standard QM/MM methods, using the TPSS and B3LYP functionals, combined with the def2-SV(P) basis set. In addition, big-QM/MM calculations were performed to accurately estimate the effect of the surroundings. Several important results have been obtained. First, we show that the singlet state is appreciably more stable than the triplet state. For the minimal model 1, the energy difference is small and for structures optimised with B3LYP, the triplet state was actually the most stable. However, with the more realistic models 2 and 3, the triplet state was 26-46 kJ mol À1 less stable than the 1 Ni-H 2 state. Including the entire protein with the big-QM/MM approach, the difference increased to 78 kJ mol À1 , which is so large that the triplet state can be discarded. Thus, the triplet state is even less stable for the H 2 adduct than for the Ni-SI a state, for which our previous studies with CCSD(T), CASPT2, and RASPT2 methods indicated that the singlet state was 24 kJ mol À1 more stable than the triplet state in the protein. 25 Second, we show that H 2 prefers to bind to Ni, rather than to Fe. For the singlet state, we find both binding modes in QM/MM optimisations with both small (30-atom) and large (298-atom) QM regions, whereas for the triplet state, only the Fe-H 2 state is found. For the isolated models 1-3, the 1 Fe-H 2 state is 9-39 kJ mol À1 less stable than the 1 Ni-H 2 state. In the protein, the difference increases to 47 kJ mol À1 . This finally settles the old controversy of the binding mode of H 2 to [NiFe] hydrogenase (cf. Table 1). Our conclusion is supported by the crystal structure of CO inhibited [NiFe] hydrogenase and the suggested H 2 binding paths. [15][16][17] Third, we have compared the performance of two DFT methods, TPSS and B3LYP, with the CCSD(T) and DMRG-CASPT2 results. For the energy difference between the two singlet binding modes ( 1 Fe-H 2 and 1 Ni-H 2 ), both functionals give accurate results with errors of less than 12 kJ mol À1 (B3LYP slightly better, with a mean absolute deviation, MAD, of 4 kJ mol À1 , compared to 7 kJ mol À1 for TPSS). This is actually similar to the accuracy of DMRG-CASPT2. However, for the singlet-triplet energy splitting, the errors are much larger. TPSS is appreciably more accurate than B3LYP with MAD and maximum errors of 10 and 21 kJ mol À1 , respectively, compared to 58 and 69 kJ mol À1 for B3LYP. Thus, we recommend the TPSS method for [NiFe] hydrogenase energetics. DMRG-CASPT2 reproduce the CCSD(T) energies even better with MAD and maximum errors of 7 and 13 kJ mol À1 .
With all these important points settled regarding the H 2 complex of [NiFe] hydrogenase, we are now ready to study the full reaction mechanism of the enzyme with QM/MM methods.