Accurate calculations of geometries and singlet-triplet energy differences for active-site models of [NiFe] hydrogenase.

We have studied the geometry and singlet-triplet energy difference of two mono-nuclear Ni(2+) models related to the active site in [NiFe] hydrogenase. Multiconfigurational second-order perturbation theory based on a complete active-space wavefunction with an active space of 12 electrons in 12 orbitals, CASPT2(12,12), reproduces experimental bond lengths to within 1 pm. Calculated singlet-triplet energy differences agree with those obtained from coupled-cluster calculations with single, double and (perturbatively treated) triple excitations (CCSD(T)) to within 12 kJ mol(-1). For a bimetallic model of the active site of [NiFe] hydrogenase, the CASPT2(12,12) results were compared with the results obtained with an extended active space of 22 electrons in 22 orbitals. This is so large that we need to use restricted active-space theory (RASPT2). The calculations predict that the singlet state is 48-57 kJ mol(-1) more stable than the triplet state for this model of the Ni-SIa state. However, in the [NiFe] hydrogenase protein, the structure around the Ni ion is far from the square-planar structure preferred by the singlet state. This destabilises the singlet state so that it is only ∼24 kJ mol(-1) more stable than the triplet state. Finally, we have studied how various density functional theory methods compare to the experimental, CCSD(T), CASPT2, and RASPT2 results. Semi-local functionals predict the best singlet-triplet energy differences, with BP86, TPSS, and PBE giving mean unsigned errors of 12-13 kJ mol(-1) (maximum errors of 25-31 kJ mol(-1)) compared to CCSD(T). For bond lengths, several methods give good results, e.g. TPSS, BP86, and M06, with mean unsigned errors of 2 pm for the bond lengths if relativistic effects are considered.


Introduction
The hydrogenases are a group of enzymes that catalyse the reversible conversion of H 2 to protons and electrons. In nature, three types of hydrogenases are known, depending on the constitution of the active site, [Fe], [FeFe], and [NiFe] hydrogenases. 1 The active site of the [NiFe] hydrogenases contains a Ni ion that is coordinated to four cysteine side chains 1,2 as can be seen in Fig. 1. Two of these bridge to the Fe ion, which also has one CO and two CN À ligands.
The [NiFe] hydrogenases have been extensively studied by both experimental and theoretical methods. [1][2][3][4][5][6] Several distinct states of the enzyme have been characterised, depending on the oxidation state of the metal ions and the coordination of Fig. 1 The active site of [NiFe] hydrogenase. 65 The Ni ion is to the left, the Fe ion to the right. extraneous ligands to the bimetal site. 1 The current consensus is that the Fe ion remains in the low-spin Fe 2+ state throughout the catalytic cycle, whereas the Ni ion cycles between the Ni 2+ and Ni 3+ states, although the Ni + state may also be involved.
However, the spin state of the Ni 2+ ion is not clear. In small inorganic complexes, low-spin Ni 2+ (which is d 8 ) normally forms square-planar complexes owing to the Jahn-Teller effect. 7,8 However, with bulky ligands, high-spin tetrahedral complexes can be obtained. 9 Crystal structures of [NiFe] hydrogenases show a Ni structure that is intermediate between square planar and tetrahedral with S-Ni-S angles of 701-1701 and a twist angle between the planes formed by the Ni ion and the two bridging cysteine sulphur atoms or the Ni ion and the two terminal cysteine sulphur atoms (f) of B701 (f = 01 for a perfect square-planar structure and 901 for a tetrahedral structure). 2,10 Several experimental studies have been performed, but they give varying results for the spin state of the Ni ion: magnetic circular dichroism (MCD) spectroscopy, 11 electron paramagnetic resonance (EPR) spectroscopy 12 and magnetic saturation techniques 13 suggest a low-spin state, whereas L-edge X-ray absorption spectroscopy 14 indicates a high-spin state.
The situation is similar in theoretical calculations: various flavours of density functional theory (DFT) give different ground states. For example, hybrid functionals like B3LYP suggest a high-spin state, whereas pure functionals like BP86 favour the low-spin state. 10,[15][16][17][18][19] This is a major obstacle for theoretical studies of the [NiFe] hydrogenases and there is an urgent need to assess the quality of the various functionals. Multiconfigurational methods are promising candidates for such an application. 20,21 A few years ago, Jayapal et al. compared the BP86 and B3LYP DFT functionals to multireference Møller-Plesset second-order perturbation theory (MRMP2 22 ) calculations of the Ni-SI a state of [NiFe] hydrogenase and concluded that BP86 is the best functional to describe such systems. 18 Unfortunately, they used only small active spaces (four or six active orbitals) for their MRMP2 calculation, which may make the results unreliable.
To check this possibility, we present in this paper multiconfigurational calculations on three different models of the [NiFe] hydrogenase with appropriate active spaces, according to the experience gained from other transition-metal complexes. 23, 24 We also study how well different DFT functionals reproduce the obtained geometries and the singlet-triplet energy difference. When available, experimental data, as well as accurate coupled-cluster calculations with single, double and (perturbatively treated) triple excitations (CCSD(T)), are also used to confirm the quality of the calculations.

Methods
Three small models were studied, which are shown in Fig. 2 4 ] 2À was studied in C 2 symmetry for both the singlet and triplet states. The [Ni(edt) 2 ] 2À complex was studied in D 2 and C 2 symmetry for the singlet and triplet geometries, respectively, whereas the NiFe model was studied in C s symmetry for both states.
Calculations were performed with second-order multiconfigurational perturbation theory on a complete activespace wavefunction (CASPT2) or on a restricted active-space wavefunction (RASPT2). ANO-RCC basis sets 26 of three different sizes were employed, as described in Table 1. BS1 was used for the geometry optimizations, whereas all three basis sets were used for energy calculations. Scalar relativistic effects were included using the Douglas-Kroll-Hess approach to second order. 27,28 Cholesky-derived auxiliary basis sets acCD 29 with a default threshold of 10 À4 a.u. was used for the geometries whereas the full Cholesky approach with a threshold of 10 À6 a.u. was used for energies, to speed-up the calculations without sacrificing the accuracy. 30,31 A parallel implementation of the CASPT2/RASPT2 method was used to further speed-up the largest calculations. 32 For the [Ni(SH) 4 ] 2À and [Ni(edt) 2 ] 2À models, the active space consisted of 12 electrons in 12 orbitals, CAS (12,12). This active space contains all Ni 3d orbitals and a second set of correlating d orbitals (which is often referred to as 3d 0 ) to account for the so-called double-shell effect. 33 In addition, two bonding orbitals between Ni and the ligands were included. In the triplet state, there are two bonding orbitals (corresponding to the two singly occupied Ni 3d orbitals) that may be involved in covalent  metal-ligand bonding (in particular in the tetrahedral structure), whereas in the singlet structure there is only one such orbital, i.e. the one corresponding to the unoccupied Ni 3d orbital. Including an extra bonding orbital in the calculation of the singlet state is not straightforward, because this orbital tends to rotate out of the active space in favour of some other orbital, e.g. a core orbital. This problem was solved by first calculating the triplet state (at the singlet structure) with CAS (12,12), and then keeping the concerned bonding orbital fixed (with the supsym option) in the active space when calculating the singlet state. It should be noted that the CASPT2(12,12) energy obtained in this way for the singlet state matches the CASPT2 energy obtained for the same state with a CAS(10,11) reference wavefunction (i.e. with the extra bonding orbitals kept inactive) to within a few tenths of a kJ mol À1 . This proves that the extra bonding active orbital is completely unimportant for the singlet state. However, this additional orbital has a significant effect on the CASPT2 energy of the triplet state, affecting the singlet-triplet energy difference by up to 10 kJ mol À1 .
For the larger NiFe model, the CASPT2(12,12) calculations were compared to calculations with a larger active space, in which ten orbitals arising from the Fe centre were added, viz. the five Fe 3d orbitals, three Fe 3d 0 orbitals (i.e. only for the occupied 3d orbitals) and two bonding ligand orbitals corresponding to the unoccupied Fe 3d orbitals. This gives a total of 22 electrons in 22 orbitals, which is intractable with current CASSCF/CASPT2 technology. Therefore, RASSCF/RASPT2 was used instead. For the singlet state, the second bonding orbital on Ni was left out of the active space (making use of the experience from the calculations on the two smaller models, showing that including the second ligand orbital does not alter the CASPT2 energy of the singlet state to any significant extent), thus giving a global RAS (20,21) space. The global active space of the RASSCF calculations was further subdivided as RAS (20,4,4;9,2,10) for the singlet versus RAS (22,4,4;8,4,10) for the triplet state, i.e. with the empty or singly occupied Ni 3d orbitals as well as their bonding ligand counterpart in the RAS2 space. In these RASSCF calculations, up to quadruple excitations were allowed out of RAS1 and into RAS3. This is essential to describe the combined double-shell effect on both metal centres. 34,35 A plot of the active orbitals obtained from the RAS(20,4,4;9,2,10) calculation on the NiFe singlet state is provided in Fig. 3. Similar plots for the other models and states are included in Fig. S1-S3 in the ESI. † For all three models, CCSD(T) 36 calculations were also performed on structures optimised with the CASPT2 ([Ni(SH) 4 ] 2À and [Ni(edt) 2 ] 2À ) or M06 (NiFe model) methods, using the ANO-RCC basis sets BS1 and BS2.
Moreover, the results were compared to calculations with DFT methods. Four DFT methods were studied with the BS1 basis set: B3LYP 37 and BP86, 38,39 as in the studies by Jayapal and Bruschi, 17,18 as well as the M06 and M06-L functionals, 40 which have been specifically parametrised for transition metal complexes.
In the CASPT2, RASPT2, and CCSD(T) calculations, all valence electrons were correlated, including the Ni and Fe 3s and 3p electrons. The CASPT2 calculations were performed with the standard IPEA (0.25-)shifted zeroth-order Hamiltonian 41 and an imaginary level shift of 0.1 a.u. The CCSD(T) calculations on the triplet state were performed starting from a ROHF wavefunction, and with restrictions on the amplitudes, such that the linear part of the wavefunction becomes a spin eigenfunction. 42 The CASPT2 and RASPT2 calculations, as well as the DFT calculations with basis set BS1, were performed using the MOLCAS 43 program  package, whereas for the CCSD(T) calculations we used the MOLPRO software. 44 For the [Ni(SH) 4 ] 2À and [Ni(edt) 2 ] 2À models, full geometry optimizations were performed with all methods (except CCSD(T), which employed the CASPT2 structures) and the BS1 basis set. The CASPT2 optimisation employed numerical gradients. However, a geometry optimization of the NiFe model at the RASPT2 level was too costly. Therefore, based on the results for the other two models, the RASPT2 and CCSD(T) calculations were performed on a structure optimised with the M06 functional.
In addition, ten DFT functionals were tested using the Turbomole software 45 50 TPSS, 51 B97-D, 52 TPSSH, 53 B3LYP, 37 PBE0, 54 and BHLYP. 55 The def2-TZVP basis set has been shown to give singlet-triplet energy differences within 1 kJ mol À1 of those obtained with the larger cc-pVQZ basis set. 10 The calculations with the semi-local functionals were sped up by expansion of the Coulomb interactions in auxiliary basis sets, the resolution-of-identity (RI) approximation. 56,57 For three of the methods (BP86-D, B3LYP-D, and B97) we also tested to include the DFT-D2 empirical dispersion correction in the calculations (both for geometries and energies). 58,59 Finally, we also studied some models with relevance to the protein structure. First, we employed singlet and triplet structures of the Ni-SI a state in [NiFe] hydrogenase from Desulfovibrio fructosovorans, obtained with quantum refinement (X-ray crystallographic refinement supplemented by combined QM and molecular mechanics, QM/MM, calculations 60,61 ). 62 From these, we cut out either the normal NiFe model, or the larger [(SCH 3 ) 2 Ni(SCH 3 ) 2 Fe(CO)(CN) 2 ] 2À model. Second, we calculated standard QM/MM structures (without restraints to the X-ray structure) for the singlet and triplet structures of the Ni-SI a , using the ComQum software. 63,64 They were based on the same crystal structure of [NiFe] hydrogenase from Desulfovibrio fructosovorans, 65 after addition of protons and solvation water molecules, and after a molecular dynamics equilibration of the whole structure. 66 The QM/MM calculations employed the def2-SV(P) basis set, 46 the Amber 99SB force field, 67 and either the TPSS or BP86 functionals. In these calculations, the QM system consisted of the [(SCH 3 ) 2 Ni(SCH 3 ) 2 Fe(CO)(CN) 2 ] 2À model, enhanced by a protonated (neutral) acetate group, as a model of Glu-25, which forms a hydrogen bond to the Cys-543 ligand, and a protonated (positively charged) imidazole group as a model of His-79, which forms a hydrogen bond to the Cys-546 ligand, in total 46 atoms. 62 In one set of TPSS calculations, this QM system was enlarged to 763 atoms (single-point energy) by including all atoms within 4.5 Å of the central [(SCH 3 ) 2 Ni(SCH 3 ) 2 Fe(CO)(CN) 2 ] 2À model, all buried charged groups connected to the active site, and moving any cut bonds at least three residues away from the active site. 66,68 3 Results and discussion In this paper, we study three simple Ni complexes (shown in Fig. 2) related to the active site in [NiFe] hydrogenase with three accurate wavefunction methods, CASPT2, RASPT2, and CCSD(T). We also compare the results with those obtained with DFT methods. We start by discussing geometries, especially those of the [Ni(edt) 2 ] 2À complex, for which accurate experimental data are available for the singlet state. 25 Then, we compare the singlet-triplet energy difference, which is expected to be even more sensitive to the level of theory. 10,17

Geometries
The optimised structures of the singlet state of the [Ni(edt) 2 ] 2À complex have D 2 symmetry (even when started from an asymmetric structure and optimised without symmetry constraints). Therefore, there is only one distance of each kind. On the other hand, the experimental structure 25 has two slightly different Ni-S bonds (219.1 and 219.8 pm). The uncertainties in the experimental distances are 0.1, 0.6, and 0.8 pm for the Ni-S, S-C, and C-C bond lengths respectively. Fig. 4 shows how well the optimised structures reproduce the experimental data. It can be seen that the CASPT2 method reproduces the experimental geometry excellently, with errors of only B1 pm for all three distances. Among the DFT functionals, all methods except SVWN overestimate the bond lengths, by up to 9 pm. The BP86 and B3LYP calculations performed using the MOLCAS software and the BS1 basis set give slightly better results than the corresponding Turbomole calculations. This is caused by differences in the basis set: turning off the Douglas-Kroll-Hess treatment of the relativistic effects changes the bond lengths by only 0.3 pm (for B3LYP), whereas shifting to the ANO-L basis set of equal size changes the three distances by 0.3-2.8 pm, increasing the mean unsigned error (MUE) from the experimental structure to 4.5 pm, i.e. close to the results obtained with Turbomole. However, the only difference between the ANO-L and ANO-RCC basis sets is that the latter was designed for relativity, so in the end, the difference is caused by relativistic effects.
The DFT-D2 dispersion correction gives very small changes (up to 0.3 pm), which do not improve the results. The SVWN and M06 methods give results closest to experiments, with MUEs of 1.7 pm. BLYP gives the largest MUE of 6 pm and BHLYP gives the largest individual difference of 9 pm for the Ni-S distance.
For the triplet state, no experimental results are available. Based on the results for the singlet, we use the CASPT2 results as the reference instead. The optimised triplet state has C 2 symmetry, so there are two distinct instances of each distance. For the triplet, the DFT methods give more varying results, with errors in the Ni-S distances from À10 pm for SVWN to 10 pm for BHLYP (Fig. 5). TPSS gives the smallest errors, less than 1 pm for both Ni-S distances. For the S-C and C-C distances, the variation is much smaller, up to 3 pm. However, for the C-C distance, all Turbomole calculations, except that with SVWN, give too long bonds of 2-3 pm. On average, the M06-L method gives the smallest errors with a MUE of 1 pm. The SVWN and BHLYP methods give the largest MUE of 4 pm.
Finally, we note that the DFT results for both spin states are somewhat different from those found by Bruschi et al. 17 (our Ni-S bonds are B4 pm shorter). The reason for this is that we use slightly different basis sets (although both are of triplezeta quality).
We have also optimised the structure of the [Ni(SH) 4 ] 2À model with the CASPT2 method and the various DFT methods in the singlet and triplet states. Both structures have C 2 symmetry, so there are two distinct Ni-S and two S-H distances, although for the singlet, the difference is less than 0.02 pm. The results are presented in Table S1 in the ESI. † As for the [Ni(edt) 2 ] 2À model, SVWN gives the best structure for the singlet state with a MUE of 2 pm compared to the CASPT2 structure. BLYP and BHLYP give the largest MUEs of 6 pm and a maximum error of 11 pm for the Ni-S distance. For the triplet state, M06 gives bond lengths that all agree with CASPT2 within 1 pm. DFT calculations using MOLCAS give significantly smaller errors than the Turbomole calculations, which all give errors of 5-12 pm for at least one Ni-S distance. As usual, SVWN and BHLYP give the largest errors.
Taking all four complexes together, BP86 calculated using MOLCAS gives the smallest error in the Ni-S distances (MUE = 2.5 pm) and also the smallest maximum error of 4 pm. M06 and M06-L give the smallest average error for all bonds (MUE = 2 pm), as can be seen in Fig. 6. Among the Turbomole calculations, TPSS gives the best Ni-S distances (MUE = 3.2 pm) and also the smallest maximum error (5 pm). PBE0 gives the smallest average error for all bonds (MUE = 2.5 pm), but several  Bond lengths for the NiFe model (for which we do not have any reliable reference structures) are presented in Table S2 (ESI †).

Singlet-triplet energy difference
Next, we looked at the energy difference between the lowest singlet and triplet states (DE ST = E triplet À E singlet ). We calculated this energy difference both vertically (for the singlet geometry) and adiabatically (for the optimum geometries of the singlet and triplet states, respectively). Table 2 shows the results for the [Ni(SH) 4 ] 2À model obtained with CASPT2 and CCSD(T) with three different basis sets (BS1, BS2, and BS3, described in Table 1).
It can be seen that the basis-set dependence is modest. The vertical DE ST changes by only B4 kJ mol À1 and the adiabatic energies change by 11-12 kJ mol À1 when going from the triple-z BS1 to the quadruple-z BS3. In particular, the mixed BS2 basis set reproduces the BS3 results within 3 kJ mol À1 .
Our best estimates of the adiabatic DE ST for the [Ni(SH) 4 ] 2À model are 12 kJ mol À1 for CASPT2 and 7 kJ mol À1 with CCSD(T), which is a reasonable agreement. For the vertical energy, the results are somewhat more different, 90 and 99 kJ mol À1 , respectively. Most importantly, both sets of calculations indicate that the [Ni(SH) 4 ] 2À model has a singlet ground state.
The results for the [Ni(edt) 2 ] 2À model are similar, as can be seen in Table 3. The basis-set dependence is modest: when the basis set is improved, DE ST increases by 4-10 kJ mol À1 . In particular, the results obtained with BS2 and BS3 agree within 2 kJ mol À1 , which is appropriate because for this larger model, we could not perform the CCSD(T) calculation with the BS3 basis set. Both the CASPT2 and CCSD(T) methods conclusively predict a singlet ground state for this model, in agreement with experiments. 25 The two methods give DE ST energies that agree within 9-12 kJ mol À1 . As for the [Ni(SH) 4 ] 2À model, CCSD(T) gives a larger difference for the vertical energy and a smaller difference for the adiabatic energy, compared to CASPT2.   In order to check the accuracy of the RASPT2 approach (which will be used for the NiFe model), test calculations with this method were performed for [Ni(edt) 2 ] 2À with the BS2 basis set. We used a RAS2 composition of two orbitals for the singlet and four orbitals for the triplet state (cf. the Methods section) and allowed for either up to double (SD) or quadruple (SDTQ) excitations out of RAS1 and into RAS3. As can be seen from Table 3, DE ST is considerably smaller with RASPT2(SD) than with CASPT2 (by up to 17 kJ mol). With RASPT2(SDTQ) DE ST is significantly increased, but it remains smaller than CASPT2 by 6 kJ mol À1 . One may therefore expect a similar underestimation of DE ST with RASPT2 for the NiFe model.
Finally, we also studied the binuclear NiFe model. For this model, both CASPT2 and RASPT2 calculations were performed. For the CASPT2 calculations, the same active space was used as for the Ni-only models, i.e. not including any orbitals centred on Fe. In the RASPT2 calculations, a set of ten orbitals (with ten electrons) on the Fe centre was added to the active space. The subdivision of the RASPT2 active space was the same as that for the [Ni(edt) 2 ] 2À mode and is detailed in the Methods section.
The results of these calculations are listed in Table 4. As before, the basis-set dependence is small and the BS2 basis set gives results that agree with those of the BS3 basis set to within 2 kJ mol À1 . All three methods conclusively predict a singlet ground state for the NiFe model, by 61 (CASPT2), 48 (RASPT2), and 57 kJ mol À1 (CCSD(T)). Somewhat strikingly, the CASPT2 results agree more closely with CCSD(T) than the RASPT2 results. In fact, a similar agreement between CASPT2 and CCSD(T) is found as for the smaller complexes, with CASPT2 predicting a slightly smaller vertical DE ST and a slightly larger adiabatic DE ST than CCSD(T). On the other hand, the RASPT2 results are systematically lower than CASPT2 (by 11-13 kJ mol À1 ). Given that RASPT2 was found to underestimate DE ST also for the [Ni(edt) 2 ] 2À complex, the CASPT2 results in Table 4 might be considered to be more accurate than the RASPT2 results. Apparently, non-dynamical correlation effects on the Fe centre are similar for the singlet and the triplet state, so that the quality of their treatment (in the active space or in the perturbational step) does not significantly affect the relative energy of the two states. This is related to the fact that the singlet-triplet transition is essentially located on the Ni-centre.
The fact that the difference between RASPT2 and CASPT2 is larger for the NiFe model, B12 kJ mol À1 , than for the [Ni(edt) 2 ] 2À complex, 6 kJ mol À1 , might be related to the fact that with two metal centres, a higher-than-four excitation level might be necessary in RASPT2 to obtain the same degree of accuracy compared to the RASPT2(SDTQ) calculation of the Ni-only complex. Therefore, we performed a RASPT2 calculation on the NiFe complex with BS2, in which the excitation level in the RASSCF wavefunction was increased to five. DE ST obtained from this RASPT2(SDTQ5) calculation is 50 kJ mol À1 , i.e. 2 kJ mol À1 higher than from the corresponding RASPT2(SDTQ) calculation. A slight further increase might be expected from the inclusion of sextuple and higher excitations in the RASSCF reference wavefunction.

Multiconfigurational character of the models
An important question when interpreting the results is how large the multiconfigurational character is for these models. Looking at the [Ni(edt) 2 ] 2À model as an example (with the BS2 basis set), we find that at the Hartree-Fock level, the triplet state (at its optimum structure) lies below the singlet state by 266 kJ mol À1 , whereas at the correlated level their relative energies are reversed and the triplet energy becomes higher than the singlet energy by 50-60 kJ mol À1 . This means that there is a huge correlation effect on DE ST , 4300 kJ mol À1 . The fact that CCSD(T) and CASPT2 differ in their description of this correlation energy by B10 kJ mol À1 is therefore not surprising, nor is the fact that very large basis sets are needed to accurately describe the relative correlation energy of both states.
At the CASSCF (12,12) level, the triplet state still lies below the singlet state by 147 kJ mol À1 . An important part of the relative correlation energy may therefore be captured by a multiconfigurational wavefunction including just a limited number of orbitals and configurations. This would indicate that this is indeed a multiconfigurational problem. Other signs are that MP2 completely breaks down, that the effect of triples on the coupled-cluster results is substantial (at the CCSD level, DE ST = À8 kJ mol À1 ), and that the effect of spin correction on the CCSD(T) energy of the triplet state is substantial. Starting from an appropriate triplet ROHF wavefunction, the difference between the UCCSD(T) and RCCSD(T) energy of the [Ni(edt) 2 ] 2À complex (BS2) amounts to 12 kJ mol À1 .
On the other hand, it is likely that this kind of multiconfigurational problem can be properly treated by CCSD(T). For both states, the occupation numbers of the natural orbitals resulting from the CASSCF calculations remain rather close to either zero or two (within 0.09; cf. Fig. 3 and Fig. S1-S3, ESI †). This is also supported by the fact that CASPT2 and CCSD(T) give results that agree to within 3-12 kJ mol À1 in all cases, which is reasonable given the huge correlation energy difference that needs to be captured.
Criteria to judge the multiconfigurational character, and hence the accuracy to be expected from CCSD(T) for transition-metal-containing molecules were recently proposed by Jiang et al. 69 The following diagnostic criteria were proposed for the computation of reliable d-block energetic and spectroscopic properties using single-reference-based model chemistries:  represent the Frobenius norm and matrix 2-norm of the coupled-cluster amplitudes for single excitations, respectively, whereas |%TAE| stands for the percentage of the (T) contribution to the total atomization energy. A too high |%TAE| value indicates that the effect of the triple correction to the CC results is substantial and hence that CCSD(T) may fail. The values of T 1 , D 1 and |%TAE| obtained from the CC calculations are shown in Table 5. For the two small models, [Ni(SH) 4 ] 2À and [Ni(edt) 2 ] 2À , one might safely conclude that non-dynamical correlation effects are limited, even though the D 1 diagnostics are rather large. This indicates that the CCSD(T) results for these two complexes should be reliable. On the other hand, for the NiFe model, both the D 1 and T 1 diagnostics are above the limits although the values of |%TAE| are still below 10%. This shows that non-dynamical correlation effects become more important in the bimetal complex, putting some doubt on the CCSD(T) results. However, this is not necessarily problematic, because the extra non-dynamical correlation effects in the NiFe complex are in fact situated on the Fe centre, and have been shown (in the discussion of CASPT2 versus RASPT2 above) not to significantly affect the relative energy of the singlet-triplet transition, which is essentially localised on the Ni centre.

Comparison with DFT results
We have also calculated DE ST for the three complexes with 14 different DFT methods (always on geometries optimised with the same method). The results are presented in Fig. 7 and they are presented as differences with respect to the best CCSD(T) results.
It can be seen that SVWN always overestimates the stability of the singlet state, by 11-73 kJ mol À1 . On the other hand, the hybrid functionals overestimate the stability of the triplet state by 14-165 kJ mol À1 . The error is roughly proportional to the amount of Hartree-Fock exchange in the functional (10% for TPSSH, 20% for B3LYP, 25% for PBE0, and 50% for BHLYP). The only exception is M06, which contains 27% exact exchange but gives a MUE of only 25 kJ mol À1 and actually gives a too high adiabatic DE ST for the [Ni(SH) 4 ] 2À model. However, all the semi-local functionals (except M06-L) give lower errors, with MUEs of 12-18 kJ mol À1 . All these methods give varying signs of the errors for the six energies, although they always overestimate the stability of the triplet state in the singlet geometry. The best results are obtained for the BP86, TPSS, and PBE functionals, which have MUEs of 12-13 kJ mol À1 . TPSS gives the smallest MUE for the adiabatic energies (7 kJ mol À1 ), whereas PBE gives the smallest MUE for the vertical energies (19 kJ mol À1 ) and also the smallest maximum error (25 kJ mol À1 ; the other semi-local functionals have maximum errors of 26-38 kJ mol À1 ).  ; the actual errors are À136, À139, À140, À161, À165, À149, and 148 kJ mol À1 for the seven entries in the legend, respectively.

Paper PCCP
The dispersion correction has a minimal effect (up to 3 kJ mol À1 ) for both BP86 and B3LYP. The MOLCAS calculations with the larger BS1 basis set give results similar to those of the Turbomole calculations, with differences of up to 9 kJ mol À1 . The MUEs agree within 1 kJ mol À1 .
Our results agree with those of the previous study by Bruschi et al. 17 in that the BP86 functional favours the singlet state and that B3LYP favours the triplet state for the adiabatic energies. However, even with B3LYP, we find a singlet ground state for all three model complexes, in contrast to the previous prediction that the triplet state should be 4 kJ mol À1 more stable than the singlet state for the [Ni(edt) 2 ] 2À model. In our study, only the PBE0 and BHLYP functionals give a triplet ground state for this complex (in disagreement with the experimental data 25 ). Only the BHLYP method gives a triplet ground state for the NiFe complex, whereas all hybrid functionals except M06 predict a triplet ground state for [Ni(SH) 4 ] 2À .

The spin state in the protein
We have also performed calculations on the NiFe model in the geometry observed inside the protein, employing structures from quantum-refinement calculations (X-ray crystallographic refinement supplemented by QM/MM calculations 60,61 ) for [NiFe] hydrogenase from Desulfovibrio fructosovorans. 62 Owing to the restraints from the crystallographic raw data, both the singlet and triplet structures are far from square planar around the Ni ion, as can be seen in Fig. 8 (the f twist angle is 68 and 711 in the two states, respectively). Consequently, the singlet state is strongly destabilised with respect to the triplet state in these structures. Table 6 shows DE ST for these structures, calculated with CASPT2, RASPT2, and CCSD(T). It can be seen that in all cases, the triplet state is the most stable by 5 (CCSD(T)) to 28 (RASPT2) kJ mol À1 . TPSS gives a result close to that of CCSD(T), 9 kJ mol À1 in favour of the triplet state.
Next, we also tested the slightly larger (and more realistic) [(SCH 3 ) 2 Ni(SCH 3 ) 2 Fe(CO)(CN) 2 ] 2À model, with methyl groups on the four sulphur ligands. The results for these calculations (only CASPT2 and RASPT2) are also included in Table 6 and show that these methyl groups stabilise the singlet state by B20 kJ mol À1 . Consequently, CASPT2 predicts that the two states are degenerate and by extrapolation towards the CCSD(T) results (which could not be calculated), we would predict that the complex in the protein should have a singlet ground state by B15 kJ mol À1 . TPSS calculations give DE ST = 4 kJ mol À1 .
These structures were obtained with a crystal structure of an oxidised form of the enzyme. 65 Unfortunately, most hydrogenase structures have problems with oxidations and photoreduction, giving structures that are a mixture of various spectrometric states. 62 Therefore, we have also studied structures optimised by combined quantum mechanics and molecular mechanics (QM/MM) for the Ni-SI a state. 66 These structures are still far from square planar, but they are less tetrahedral around the Ni ion with f twist angles of 58 and 651-681, respectively. However, they do not change the relative stability of the two states very much. At the TPSS level and using a 46-atom QM system (with methyl groups on the sulphur ligands and including a protonated acetate group and a protonated imidazole group, forming hydrogen bonds with two of the sulphur ligands), DE ST = 9 kJ mol À1 for structures optimised with the BP86 functional, but À3 kJ mol À1 for structures optimised with the TPSS functional (Table 7). If a point-charge model of the surrounding protein is included in the QM calculation, the singlet state is stabilised by 3-7 kJ mol À1 ,  but if also the MM terms are added, yielding a full QM/MM energy, DE ST = 5 or À4 kJ mol À1 , respectively. All these energies were obtained with the smaller def2-SV(P) basis set. If we instead use the larger def2-TZVP basis set (used in the other DFT calculations in this article), the singlet state is stabilised by 9-14 kJ mol À1 . Finally, we calculated DE ST for a very large 763-atom QM system, including all atoms within 4.5 Å of the central [(SCH 3 ) 2 Ni(SCH 3 ) 2 Fe(CO)(CN) 2 ] 2À model, all buried charged groups connected to the active site, and moving any cut bonds at least three residues away from the active site, and including a point-charge model of the surrounding protein and solvent. 66,68 This gave DE ST = 6-11 kJ mol À1 , which should include all important effects from the QM/MM calculation. Therefore, we arrive at our final estimate DE ST = 24 kJ mol À1 for both structures by adding the correction for going from the def2-SV(P) to the def2-TZVP basis set and the difference between TPSS and CCSD(T) for the NiFe model in the crystallographic structure. This indicates that the active site is most likely in the singlet state in the protein, but the energy difference is so small that the conclusion is still somewhat uncertain.
Jayapal et al. also studied the [(SCH 3 ) 2 Ni(SCH 3 ) 2 Fe(CO)(CN) 2 ] 2À model with the MRMP2 method and found that the singlet state was 60-99 kJ mol À1 more stable than the triplet state for structures optimised in vacuum, 18 i.e. similar to our results, considering that their model included methyl groups on the thiolate ligands. They used a much smaller active space than in our calculations with only four orbitals: a pair of bonding and antibonding orbitals between a Ni 3d orbital and the four surrounding sulphur atoms, an essentially pure Ni 3d orbital, and a diffuse orbital that is the correlating 3d 0 orbital of the latter (they also performed a test calculation with six active orbitals, selected according to the occupation numbers from an unrestricted Kohn-Sham calculation). CASPT2 calculations with such a (4,4) active space on the NiFe model (with BS2) gave vertical and adiabatic DE ST of 98 and 43 kJ mol À1 , respectively, i.e. 17-19 kJ mol À1 lower than the corresponding CASPT2 (12,12) calculations shown in Table 4, and 31 or 14 kJ mol À1 lower than the corresponding CCSD(T) results. This shows that the current results are more accurate.

Conclusions
In this article, we have studied the accuracy of the CASPT2 and RASPT2 methods, as well as of 14 DFT functionals, for three models of the [NiFe] hydrogenase active site (shown in Fig. 2). The CASPT2 and RASPT2 approaches with appropriate active spaces and large basis sets accurately reproduce experimental and CCSD(T) results (bond lengths within 1 pm, singlet-triplet energy differences within 15 kJ mol À1 ). Therefore, the CASPT2 and RASPT2 methods can be used as a suitable reference when there are no experimental data and when CCSD(T) calculations become intractable or unreliable because of strong non-dynamical correlation.
Concerning the DFT methods, we have confirmed previous observations that B3LYP gives too long metal-ligand bonds. 17,70 Moreover, our results illustrate the usual tendency of semilocal functionals to over-stabilise the low-spin state and of hybrid functionals to over-stabilise the high-spin state. 17,71 The BP86, TPSS, and PBE methods give singlet-triplet energy differences in best agreement with the CCSD(T) results (MUEs of 12-13 kJ mol À1 and maximum errors of 19-29 kJ mol À1 ). TPSS also gives the most accurate Ni-S bond lengths and, together with PBE, BLYP, TPSSH, and PBE0, the best results for all bond lengths in the Turbomole calculations with the def2-TZVP basis set, so it seems to be the method of choice for this type of systems. The geometries from the DFT calculations seem to improve if relativistic effects are considered, but DE ST is not significantly changed.
Our results suggest a singlet ground state of Ni-SI a in [NiFe] hydrogenase, in agreement with previous BP86 studies, 17,18 but in contrast to B3LYP results. 15,16 For vacuum-optimised structures, DE ST is so large (57 kJ mol À1 for the NiFe model and B20 kJ mol À1 larger for the [(SCH 3 ) 2 Ni(SCH 3 ) 2 Fe(CO)(CN) 2 ] 2À model) and the methods so advanced that we strongly believe that this result can be trusted. However, the geometry observed in the protein strongly disfavours the singlet state and with this geometry, the singlet state is predicted to be only B24 kJ mol À1 more stable than the triplet state, which is on the limit of the accuracy of the methods used. Apparently, it seems that the protein has selected a structure of the active site for which the singlet and triplet states have similar energies, as has also been suggested in a recent study, showing that the minimal energy crossing point between the singlet and triplet energy surfaces is close to the structure observed in the protein. 10 In future publications, we will investigate how this selection is obtained and what mechanistic advantages it gives.