Extremely large differences in DFT energies for nitrogenase models

Nitrogenase is the only enzyme that can cleave the triple bond in N 2 , making nitrogen avaiable for other organisms. It contains a complicated MoFe 7 S 9 C(homocitrate) cluster in its active site. Many computational studies with density-functional theory (DFT) of the nitrogenase enzyme have been presented, but they do not show any consensus – they do not even agree where the first four protons should be added, forming the central intermediate E 4 . We show that the prime reason for this is that different DFT methods give relative energies that differ by almost 600 kJ mol (cid:2) 1 for different protonation states. This is 4–30 times more than what is observed for other systems. The reason for this is that in some structures, the hydro-gens bind to sulfide or carbide ions as protons, whereas in other structures they bind to the metals as hydride ions, changing the oxidation state of the metals, as well as the Fe–C, Fe–S and Fe–Fe distances. The energies correlate with the amount of Hartree–Fock exchange in the method, indicating a variation in the amount of static correlation in the structures. It is currently unclear which DFT method gives the best results for nitrogenase. We show that non-hybrid DFT functionals and TPSSh give the most accurate structures of the resting active site, whereas B3LYP and PBE0 give the best H 2 dissociation energies. However, no DFT method indicates that a structure of E 4 with two bridging hydride ions is lowest in energy, as spectroscopic experiments indicate.


Introduction
2][3] The nitrogenase reaction is essential to the life on earth -although 78% of the atmosphere is N 2 , nitrogen is often the limiting factor for plant growth and a main component of synthetic fertilisers. 3In fact, the industrial Haber-Bosch process to form ammonia consumes B1% of the world's energy supplies, produces half of the total biologically available nitrogen on earth 4 and is a major factor in the recent human population explosion. 36][7][8][9] The catalytic centre is the MoFe 7 S 9 C(homocitrate) FeMo cluster bound to the protein by a cysteine and a histidine residue.In some enzymes, the Mo ion is replaced by vanadium or iron. 10The protein also contains a Fe 8 S 7 Cys 6 cluster, called the P cluster, which transfers electrons.The electrons are provided by another protein, called the Fe protein, which binds two ATP molecules.Hydrolysis of the ATP molecules triggers the dissociation of the Fe protein, opening up for additional electron transfers.
Nitrogenase catalyses the N 2 + 8e À + 8H + + 16ATP -2NH 3 + H 2 + 16ADP + 16P i reaction.It is typically described by the Lowe-Thorneley scheme, involving nine intermediates, E 0 -E 8 , differing in the number of electrons and protons delivered to the enzyme. 112,13 Of particular interest is E 4 , which is believed to be the species that binds N 2 concomitant with the release of H 2 by reductive elimination. 3][15][16][17][18][19][20][21][22][23][24][25][26] Unfortunately, they have led to strongly diverging mechanistic suggestions.For example, some studies have proposed that N 2 is sequentially protonated first on one N atom, which dissociates into NH 3 before the second atom is protonated, 23,25 whereas other studies have suggested that it is alternatively protonated on the two N atoms, forming N 2 H 2 and N 2 H 4 intermediates. 22,26Likewise, there is no agreement on how N 2 binds to the cluster, e.g.side-on to one Fe ion, 22 with one N atom bridging two Fe ions, after the dissociation of one of the sulfide ions, 23 in the centre of the cluster, displacing a triply protonated carbide ion, 24 or forming a covalent bond to the carbide ion. 25,26n fact, there is not even any consensus regarding the structure of the central intermediate E 4 .For example, Hoffman and coworkers have proposed models with two protonated sulfide ions and two bridging hydride ions, 27 in agreement with spectroscopic results, whereas Siegbahn has argued that it is energetically much more favourable to triply protonate the central carbide ion. 28Recently, we have performed a systematic combined quantum mechanical and molecular mechanical (QM/MM) study of the possible protonation states of the FeMo cluster, providing a great wealth of possible structures for the various intermediates in the reaction mechanism. 29Here, we examine ten of these and show extreme energy differences between different DFT methods, partly explaining the discrepancy found in previous studies.We also discuss which DFT method gives the more accurate results, by studying the geometry of the resting E 0 state and H 2 dissociation from E 4 and E 2 .

The protein
All calculations are based on the 1.0 Å crystal structure of nitrogenase from Azotobacter vinelandii (PDB code 3U7Q). 70][31] The entire heterotetramer was included in the calculations and the DFT calculations were concentrated on the FeMo cluster in the C subunit.The P clusters and the FeMo cluster in subunit A were modelled by MM in the fully reduced and resting states, respectively. 30he protonation states of all residues were the same as before: 30 all Arg, Lys, Asp, and Glu residues were charged, except Glu-153, 440, and 231D (a letter ''D'' after the residue number indicates that it belongs to that subunit; if no letter is given, it belongs to subunit C; subunits A and B are identical to the C and D subunits).Cys residues coordinating to Fe ions were deprotonated.His-274, 451, 297D, 359D and 519D were assumed to be protonated on the ND1 atom, His-31, 196, 285, 383, 90D, 185D, 363D and 457D were protonated on both the ND1 and NE2 atoms (and therefore positively charged), whereas the remaining 14 His residues were modelled with a proton on the NE2 atom.The homocitrate was modelled in the singly protonated state with a proton shared between the hydroxyl group (which coordinates to Mo) and the O1 carboxylate atom.This protonation state was found to be the most stable in a recent extensive QM/MM, molecular dynamics and quantumrefinement study 30 and it is also supported by another study. 32he protein was solvated in a sphere with a radius of 65 Å around the geometrical centre of the protein.160 Cl À and 182 Na + ions were added to neutralise the protein and give an ionic strength of 0.2 M. 33 The final system contained 133 915 atoms.The added protons, counter ions and water molecules were optimised by a simulated annealing calculation (up to 370 K), followed by a minimisation, keeping the other atoms fixed at the crystal-structure positions. 30ll MM calculations were performed with the Amber software. 34For the protein, we used the Amber ff14SB force field 35 and water molecules were described by the TIP3P model. 36For the metal sites, the MM parameters were the same as in our previous investigation. 30The metal sites were treated by a nonbonded model 37 and charges were obtained with the restrained electrostatic potential method, obtained at the TPSS/def2-SV(P) level of theory 38,39 and sampled with the Merz-Kollman scheme. 40

QM/MM calculations
The QM/MM calculations were performed with the ComQum software. 41,42In this approach, the protein and solvent are split into two subsystems: system 1 (the QM region) was relaxed by QM methods, whereas system 2 contained the remaining part of the protein and the solvent.It was kept fixed at the original coordinates (equilibrated crystal structure). 30n the QM calculations, system 1 was represented by a wavefunction, whereas all the other atoms were represented by an array of partial point charges, one for each atom, taken from the MM setup.When there is a bond between systems 1 and 2, 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 atoms, CL) in the full system. 41,43All MM atoms were included in the point-charge model, except the CL atoms. 44he total QM/MM energy in ComQum was calculated as 41,42 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 system 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 region set to zero (to avoid double-counting of the electrostatic interactions).Thus, ComQum employs a subtractive scheme with electrostatic embedding and van der Waals link-atom corrections. 45he geometry optimisations were 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.
It should be noted that the M06-2X and M06-HF functionals were constructed for main-group thermochemistry and chargetransfer excitations, respectively, and are not recommended for systems with significant static correlation, like transition metals. 56The same applies also to BHLYP with its large amount of HF exchange. 63These three methods reproduce the structure of the resting state quite poorly (Fig. S3, ESI †) and they often give rise to a qualitative change in the geometry of the optimised structures of models 5-7 and 9 (Table S1, ESI †).Therefore, these three functionals were excluded from the figures in the main article, but the results are included in the ESI † (Fig. S1-S3), because they emphasize the trends with respect to % HF.
The FeMo cluster was modelled by MoFe 7 S 9 C(homocitrate) (CH 3 S)(imidazole), where the last two groups are models of Cys-275 and His-442.In addition, all groups that form hydrogen bonds to the FeMo cluster within 3.5 Å in the crystal structure 7 were also included, viz.Arg-96, His-195 and Arg-359 (sidechains), Gly-356, Gly-357 (backbone), as well as two water molecules.In total, the QM system contained 113-117 atoms (depending on the E n state, i.e. the number of added protons) and it is shown in Fig. 1.Following recent Mo ¨ssbauer, anomalous dispersion and QM investigations, 13,24,64,65 we used the oxidation stateassignment Mo III Fe II 3 Fe III 4 of the metal ions in the resting (E 0 ) state, giving a net charge of À5 for the QM system.We studied three oxidation states of the FeMo cluster, obtained by adding 0, 2 or 4 electrons to the resting states, denoted E 0 , E 2 and E 4 .One proton was added together with each electron, so the net charge of the cluster was À3 for all three states.
Experiments have shown that the ground spin state of E 0 and E 2 are quartets with a surplus of three a electrons, whereas E 4 is a doublet. 3,13These spin states were used in all calculations.The electronic structure of all QM calculations was obtained with the broken-symmetry (BS) approach: 18 for each model, we used the best of the 35 possible BS states 31 found in our previous study. 29These states are specified in Table 1, which also give further description of the various models (which atoms are protonated and the number of bonds within the cluster).The 13 calculations with different DFT methods for each complex were started with the same wavefunction, to ensure that they belong to the same BS state.Bond lengths involving the added protons are described in Table S2 (ESI †) for all models optimised with the 13 DFT methods, spin densities on the metals are Fig. 1 The FeMo cluster illustrating the QM system used in all QM/MM calculations and the naming of the atoms in the cluster.
Table 1 Description of the ten models, listing the E n state, the positions of the added protons (atom names are shown in Fig. 1; the number in brackets indicates the direction of the proton a 27 ), the BS state, b the number of bonds of various types (a Fe-Fe interaction is assumed to be broken if it is 0.5 Å longer than in the resting state; the Fe-S and Fe-C bonds are considered broken if they are 0.3 Å longer than in the resting state), and the number of Fe(II) ions in the FeMo cluster (n, formal assignment, assuming that all metal-bound hydrogen atoms are hydride ions) Undoubtedly, the present calculations involve several approximations, e.g. the use of single QM/MM structures and small basis sets without any inclusion of entropy or thermal effects.However, the extremely large difference between the DFT results will undoubtedly remain also in improved calculations.

Results
In this study, we compare energies and geometries of 10 models of the FeMo cluster in the E 4 , E 2 or E 0 state, shown in Fig. 2, calculated with 13 different DFT methods.All structures were optimised with QM/MM individually for each DFT method.The aim is to describe and understand the large variation of the results obtained with the various methods and to gain some understanding of which method gives the more reliable results.

E 4 states
We first compare the relative stability of the best protonation states of E 4 obtained with either the B3LYP (1) or TPSS (5) functionals in our extensive study, DE 15 . 29As can be seen from Fig. 2, both structures are protonated on one of the m 2 sulfide ions (S2B), whereas the other three protons are on the central carbide ion in 1 but on Fe ions in 5. Protonation of the carbide leads to significant distortions of the FeMo cluster (cf.Table 1), but the carbide ion remains inside the cluster in our QM/MM calculations, 29 in variance to the QM-cluster calculations by Siegbahn, in which the protonated carbide ion moves to the periphery of the cluster. 66rom the relative energies in Fig. 3 (blue curve), it can be seen that the (meta) GGA functionals (TPSS, PBE, BP86, BLYP, M06-L and B97D) predict that structure 5 is 8-242 kJ mol À1 more stable than structure 1.On the other hand, the hybrid functionals (TPSSh, B3LYP, PBE0 and M06) predict that structure 1 is 39-353 kJ mol À1 more stable than structure 5.The latter agrees with Siegbahn's suggestion (based on B3LYP calculations) that it is much more favourable to protonate the central carbide than the metal ions. 28However, a similar study based on a GGA functional would have reached the opposite conclusion.
Thus, the various DFT methods give DE 15 results that differ by up to 595 kJ mol À1 .8][69] For example, Fig. 4 shows the reaction energies of the N 2 + 3H 2 -2NH 3 reaction, divided into three steps (involving the formation of N 2 H 2 , N 2 H 4 and finally two NH 3 molecules).It can be seen that the 13 DFT functionals give results that agree within 35-38 kJ mol À1 .In particular, all methods agree that the first step (formation of N 2 H 2 ) is strongly uphill (by 167 kJ mol À1 on average), whereas the other two steps are downhill by 102 and 154 kJ mol À1 on average, respectively.
The hybrid functionals involve a varying amount of Hartree-Fock exchange (% HF, 10-27%), as is also shown in Fig. 3 (cyan line and right axis).There is a fair correlation between DE 15 and % HF, R 2 = 0.75, although M06 breaks the trend.In fact, we also used three functionals with a larger % HF, which increased the correlation to 0.88 and the energy range to 1097 kJ mol À1 (Fig. S1, ESI †), but since such methods are normally not recommended for transition-metal systems 56,63 and they gave qualitative changes in the structure of 5, we present these results in the ESI.† For the N 2 + 3H reactions, the correlation to the % HF is appreciably weaker, 0.03-0.42.
If the structures are not optimized with each DFT functional, the variation is even larger, 862 kJ mol À1 (using the B3LYP structure for 1 and the TPSS structure for 5; Table S3, ESI †), showing the importance of optimizing the individual structures  when the energy differences are so large.However, the correlation of DE 15 obtained with or without geometry optimization is very high, 0.98.
In 1, there are three protons on the central carbide, whereas in 5, they are instead (as hydride ions) on Fe ions.To further understand the large variation in the DFT results, we did similar calculations on three structures (2, 3 and 4) with intermediate protonation, viz.with 2, 1 and no protons on the carbide and with 0, 1 and 2 metal-bound hydride ions, respectively (cf.Fig. 2).In fact, structure 4 is one of the E 4 structures suggested by the Hoffman group (based on BP86 calculations) with two hydride ions bridging two Fe ions each and two protons on sulfide ions. 27From Fig. 3 (red line), it can be seen that DE 45 is similar for most functionals, with 5 being more stable than 4 by 5-62 kJ mol À1 .Only B3LYP and M06 give the opposite ordering (by 19-34 kJ mol À1 ).However, there is a quite strong correlation between DE 45 and % HF, R 2 = 0.8.
In model 3, there is only one hydride ion and the carbide ion is singly protonated.Interestingly, DE 34 shows an appreciably larger variation among the functionals (Fig. 3, yellow).In general, 3 is less stable than 4 according to the GGA functionals (by 35-137 kJ mol À1 ), whereas the opposite is true with the hybrid functionals (by 35-203 kJ mol À1 ).However, with B97D and M06, they are of a similar energy.There is a fair correlation between DE 34 and % HF (R 2 = 0.58; 0.76 with all 13 functionals).
Finally, model 2 has no hydride ions and a doubly protonated carbide ion.From Fig. 3 (green line), it can be seen that the GGA functionals, except M06-L and B97D indicate that 3 is 7-33 kJ mol À1 more stable than 2, whereas the other methods give the opposite result by 5-125 kJ mol À1 .There is a strong correlation between DE 23 and % HF (R 2 = 0.78 or 0.92).If we instead compare the stability of structures 2 and 1, i.e. the effect of moving a proton from a sulfide to the carbide (brown line in Fig. 3), we find that the GGA functionals and M06 indicate that structures 2 and 1 are close in energy (À10 to 2 kJ mol À1 ), whereas the other hybrid functionals indicate that structure 1 is 29-39 kJ mol À1 more stable than structure 2. There is only a low correlation between DE 12 and % HF (0.3, but it increases to 0.75 with all 13 DFT functionals).
Thus, the large variation in DE 15 among the DFT functionals is caused mainly by DE 34 and DE 23 .This indicates that the variation is related to at least two factors.The isomerizations 4 -5, 3 -4 and 2 -3 all involve a change in the number of hydride ions and therefore also the formal oxidation state of the Fe ions (a hydride ion has two more electrons than a proton; therefore, the number of Fe(II) ions in the structures 1-5 are 7, 7, 5, 3 and 1, respectively), so that structures 1 and 5 actually differ by six levels in the formal oxidation state of the metal ions.On the other hand, 3 -4, 2 -3 and 1 -2 involve the formation of C-H bonds, associated with a distortion of the FeMo cluster, leading to partial cleavage of Fe-C, Fe-S and Fe-Fe bonds.This is detailed in Table 1, in which it can be seen that eight of the Fe-Fe distances have elongated by more than 0.5 Å and four of the Fe-C bonds and one of the Fe-S bonds are elongated by more than 0.3 Å in 1, compared to 5 or the resting E 0 state.Models 2 and 3 have a lower number of broken bonds (cf.Table 1).Apparently, the variation among the DFT methods is the largest when both effects apply.

H 2 binding and dissociation
To gain some additional information about how the energies vary with the DFT functionals, we studied also the binding and dissociation of H 2 from the E 4 and E 2 states of the MoFe cluster.First, we included two complexes with H 2 bound to the FeMo cluster in the E 4 state, both with a proton on S2B.Structure 6 has H 2 bound end-on to Fe5 and the fourth H is bridging Fe2 and Fe6, whereas in structure 7, H 2 is bound side-on to Fe2 and with the fourth H on Fe4 (cf.Fig. 2; atom names are shown in Fig. 1).For all methods, end-on binding (6) was more favourable than side-on binding (7), by 20-88 kJ mol À1 (Fig. 5, blue line).Moreover, with the GGA functionals, 6 was more stable than 5 (the most stable E 4 structure with those functionals) by 7-104 kJ mol À1 (Fig. 5, red line).On the other hand, 1 was appreciably more stable than 6 for the B3LYP and PBE0 functionals (by 239-246 kJ mol À1 ), whereas the two structures were isoenergetic within 2 kJ mol À1 with TPSSh and M06.
Next, we studied two structures of the E 2 state (i.e. with two protons and electrons less), one with the two protons on the  central carbide (8) and the other with one proton on the S2B m 2 sulfide and one hydride ion bridging Fe2 and Fe6 (9; cf.Fig. 2).They are the best E 2 structures obtained in our previous study with B3LYP and TPSS, respectively. 29As expected, the GGA functionals, but also TPSSh, suggest that 9 is 3-203 kJ mol À1 more stable than 8, whereas the other hybrid functionals favour 8 by 9 (M06) or 122-153 kJ mol À1 (Fig. 5, yellow line).The correlation of DE 89 to % HF is 0.6.
Finally, we studied also the resting state of the enzyme (E 0 ; 0) to obtain an estimate of the dissociation energy of H 2 from E 2 (9/8 -0 + H 2 ).Fig. 5 (brown line) shows that it is favourable for all GGA functionals, as well as TPSSh and M06 by 88-166 kJ mol À1 , whereas it is unfavourable for B3LYP and nearly thermoneutral for PBE0.

Geometry of the resting state
For E 0 , we can also compare the QM/MM structures with the crystal structure of nitrogenase (3U7Q). 7It is at such a high resolution (1.0 Å) that it can be used to calibrate the DFT geometries.In fact, it has been used to determine the brokensymmetry state of the FeMo cluster and the protonation of the homocitrate ligand. 32Fig. 6 shows the root-mean-squared deviation (RMSD), as well as the mean absolute deviation (MAD) and maximum deviation of the metal-ligand and metal-metal distances for the various DFT methods.It can be seen that TPSSh gives the lowest RMSD (0.06 Å), whereas the best metal-ligand distances are obtained with PBE and BP86 (giving the lowest MAD = 0.02 Å) or BLYP (lowest maximum deviation, 0.06 Å).On the other hand, the best metal-metal distances are obtained with B97D and TPSSh (MAD = 0.02 Å) or BLYP methods (maximum deviation of 0.07 Å).Thus, TPSSh and the GGA functionals (except M06-L) give the best structures, appreciably better than those of the other hybrid functionals.

Spin densities
The geometry optimisations with different DFT functionals were started from the same wavefunction, to keep them as similar as possible.We used the experimentally observed spin state and the best BS state obtained for each structure in our previous investigation (listed in Table 1). 29These were obtained by B3LYP for structures 1-3 and 8, and by TPSS for the other structures.Thus, we did not attempt to find better BS states for the other DFT functionals, which would be extremely timeconsuming and would have mixed in differences caused by the use of different BS states.
The spin densities obtained for the eight metal ions in the various calculations are listed in Table S6 (ESI †).It can be seen that for most structures, the calculations with different DFT functionals remained in the same BS states, defined by the Fe ions with negative spin densities.However, for two structures, some DFT calculations gave qualitatively differences.For 6, most calculations were in the BS7-3 state, characterised with large negative spin populations on Fe3, Fe4 and Fe6 (rÀ2.4).However, the three GGA functionals PBE, BP86 and BLYP gave negative spin populations on the same three Fe ions, but that on Fe6 was small, À0.5 and Fe2 also had a negative spin of intermediate size, around À1.5 (therefore corresponding to the BS6-2 state).We could find the other state also with these functionals, but during geometry optimisation, they changed back to the original state, except for BLYP, for which it is 12 kJ mol À1 less stable than the original state.Still, this difference in the spin states is not reflected in the energies in Fig. 5.
For 4, all DFT calculations give large negative spin populations on Fe1 and Fe3.However, for five of the functionals (the GGA functionals, except M06-L), Fe7 had a small population on Fe7 (À1.2 for B97D or À0.1 for the others), whereas for the other functionals, the population is large and negative.Again, the other state could be found but it typically changed back during geometry optimisation.For BLYP, TPSSh and B3LYP, both states could be found, with energy differences of 12-97 kJ mol À1 .In this case, the difference in the wavefunction is reflected in the energies in Fig. 3, which shows an increased stability for 4 for the functionals giving a high spin population on Fe7.
In addition, there are also systematic differences between the spin densities obtained with the various methods.First, the magnitude of the spin populations increases with % HF.It is in general the lowest for BLYP, PBE and BP86 (2.5 on average) and the largest for PBE0 (3.6).The spin on the seven Fe ions also becomes more similar when % HF increases.For example, the difference between the maximum and minimum absolute spin populations of structure 1 is 1.0 for PBE and BP86, but only 0.2 for PBE0.
The difference is even larger for the spin on the Mo ion.It is negative in all structures, except 5.With the TPSS, PBE, BP86 and BLYP functionals, it is in general small in magnitude, 0.2-0.5, but increased to 0.8-1.0 for structures 1 and 2. However, for M06-L and the hybrid functionals, it is appreciably higher, by 0.6-0.8Fig. 6 Performance of the various DFT methods for the structure of the resting E 0 state compared to the crystal structure. 6For each method, the RMSD for the Mo, Fe, S and C atoms of the cluster, as well as the directly connected ligand S, O and N atoms are given, together with the MAD and maximum deviation for the 34 metal-ligand (M-L) distances and the 15 short (o3 Å) metal-metal (M-M) distances.
on average.For structures 1 and 2, these functionals indicate a high-spin state of Mo (spin populations of 1.6-2.5),[72]

Discussion
We have studied how the relative energies of various protonation and H 2 -bound states vary when they are calculated with ten different DFT methods.The results in Fig. 3 and 5 show an extremely large variation of the relative energies, up to 600 kJ mol À1 , which often correlate to % HF of the methods.A strong dependence on DFT for relative energies and a correlation to % HF have been observed previously for several types of systems.For example, Reiher et al. studied the energy difference between the high-and low-spin states for a number of Fe(II)S 4 NL complexes with L = NH 3 , N 2 H 4 , CO, NO + , PH 3 and PMe 3 . 67They showed that predictions obtained by B3LYP and BP86 differed by 94-122 kJ mol À1 and that variants of B3LYP with % HF varying between 0 and 25% gave energies that varied by up to 160 kJ mol À1 .They observed that for % HF between 10 and 15%, the results agreed with experiments and recommended the B3LYP* method with 15% HF.This is the method used by Siegbahn for energies (but not structures) for nitrogenase. 28ikewise, the homolytic Co(III)-C bond dissociation energy in methyl corrin (coenzyme B 12 without sidechains), Jensen and Ryde reported a difference of 57 kJ mol À1 between B3LYP and BP86 (51-66 kJ mol À1 for related models with different axial ligands and metals, Co, Fe or Ni). 68Later, more functionals were studied, increasing the variation to 139 kJ mol À1 (66 kJ mol À1 excluding BHLYP) and showing a correlation to % HF of R 2 = 0.89. 73or the oxy-transfer reaction in Mo-dependent dimethylsulfoxide reductase, Li et al. reported a variation of 116 kJ mol À1 for the activation energy and 103 kJ mol À1 for the reaction energy over eight DFT functionals (63 and 61 kJ mol À1 excluding BHLYP).The correlation to % HF was 0.95 and 0.98. 74The results were calibrated by using LCCSD(T0) calculations, which showed that for the activation energy, B3LYP gave the best result (9 kJ mol À1 error), whereas for the reaction energy, GGA functionals gave better results (errors of 1-6 kJ mol À1 ).Similar results were obtained if geometries were optimised with the various methods or if also W was considered. 69hus, previous studies of metalloenzymes have sometimes showed a large variation in the results of different DFT methods, especially for spin-state energies or when the oxidation state of the metal changes.However, the variation observed for nitrogenase is 4-10 times larger than what has been observed before.
Of course, the prime question is which of the DFT methods (if any) can be trusted.In general, hybrid DFT functionals give appreciably better results than GGA functionals for maingroup thermochemistry, e.g. with a weighted MAD of 21, 27 and 38 kJ mol À1 for M06-2X, B3LYP and TPSS, respectively, for the recent GMTKN55 test set (and even better results for double hybrid functionals). 75Hybrid functionals remain better for some transition-metal complexes, in particular with hard ligands. 67,76owever, there are ample examples of transition-metal complexes (especially with softer ligands), for which GGA methods give better results than hybrid functionals, e.g. for reactions related to the nitrogenase reaction, 77 for [NiFe]-hydrogenase models 78,79 and for complexes with nitride and carbonyl ligands. 80In a recent test of 41 closed-shell organometallic reactions, PBE0 gave the best results among the DFT functionals employed in the present study, but the MAD was not much higher for TPSSh and TPSS (11, compared to 12 and 14 kJ mol À1 ), whereas B3LYP was appreciably worse, 21 kJ mol À1 . 81Clearly, the performance of DFT methods depend on the system studied, probably depending on the importance of static correlation. 80n alternative is to select the functional based on comparisons with experimental observations.In this study, we have compared several properties.Spectroscopic evidence indicate that the E 4 state has two bridging hydride ions. 3Unfortunately, the results in Fig. 3 show that none of the DFT functionals suggest that the best E 4 structure is 4, which has two bridging hydride ions.
Moreover, experiments indicate that four electrons and protons should be added to the resting state before it may bind N 2 .Of course, it is then essential that H 2 does not form and dissociate before N 2 binds; otherwise two electrons and protons have been consumed with no gain.However, the energies in Fig. 5 indicate that all methods except B3LYP and PBE0 give strongly favourable H 2 dissociation energies from both E 4 and E 2 .
On the other hand, the results in Fig. 6 show that the TPSSh and GGA methods (except M06-L) best reproduce the crystal structure of the resting state.Thus, our calculations indicate that it may be hard to find a DFT method that gives good results for all properties of the FeMo cluster in nitrogenase.

Conclusions
We have shown that nitrogenase models with different protonation states give an unprecedentedly large variation in relative energies obtained with different DFT methods, up to 600 kJ mol À1 .This variation is connected to differences in the formal oxidation state of the Fe ions in the cluster (number of hydride ions), as well as to the number of Fe-C and Fe-H bonds.This explains the large discrepancy between mechanisms suggested by different groups. 27,28nfortunately, there is currently no quantum mechanical method that can be used to calibrate the DFT methods for a cluster of the size and electronic complexity as the one in nitrogenase.Instead, we suggest that energies are calculated with several DFT methods to check the sensitivity of the results.We have attempted to calibrate the methods by using experimental data.Unfortunately, GGA functionals and TPSSh give the best geometries, whereas B3LYP and PBE0 give better H 2 -dissociation energies but a high-spin Mo state for the best E 4 structures.None of the tested methods indicate that the best E 4 state has two bridging hydride groups.This calls for some humility in the interpretation of DFT results on nitrogenase.Undoubtedly, there is an urgent need of more accurate methods that can treat systems of the complexity of the FeMo cluster.

Fig. 2
Fig. 2 Structures of the 10 studied nitrogenase models, 0-9.Mo, Fe, S and C are shown in cyan, orange, yellow and grey, respectively.Added hydrogen atoms are shown in green balls.For clarity, only the central atoms in the cluster are shown.