Lili
Cao
and
Ulf
Ryde
*
Department of Theoretical Chemistry, Lund University, P. O. Box 124, SE-221 00 Lund, Sweden. E-mail: Ulf.Ryde@teokem.lu.se
First published on 21st December 2018
Nitrogenase is the only enzyme that can cleave the triple bond in N2, making nitrogen avaiable for other organisms. It contains a complicated MoFe7S9C(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 E4. We show that the prime reason for this is that different DFT methods give relative energies that differ by almost 600 kJ mol−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 hydrogens 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 H2 dissociation energies. However, no DFT method indicates that a structure of E4 with two bridging hydride ions is lowest in energy, as spectroscopic experiments indicate.
Crystallographic studies have shown that the nitrogenase is a large α2β2 heterotetramer.5–9 The catalytic centre is the MoFe7S9C(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.10 The protein also contains a Fe8S7Cys6 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 N2 + 8e− + 8H+ + 16ATP → 2NH3 + H2 + 16ADP + 16Pi reaction. It is typically described by the Lowe–Thorneley scheme, involving nine intermediates, E0–E8, differing in the number of electrons and protons delivered to the enzyme.11 Nitrogenase has been extensively studied by a great wealth of biochemical, spectroscopic and kinetic methods and several of the intermediates have been trapped and characterised.1–3,12,13 Of particular interest is E4, which is believed to be the species that binds N2 concomitant with the release of H2 by reductive elimination.3
The experimental studies have been supplemented by many density functional theory (DFT) studies, which can give an detailed atomistic and energetic picture of the reaction mechanism.3,13–26 Unfortunately, they have led to strongly diverging mechanistic suggestions. For example, some studies have proposed that N2 is sequentially protonated first on one N atom, which dissociates into NH3 before the second atom is protonated,23,25 whereas other studies have suggested that it is alternatively protonated on the two N atoms, forming N2H2 and N2H4 intermediates.22,26 Likewise, there is no agreement on how N2 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,26
In fact, there is not even any consensus regarding the structure of the central intermediate E4. 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.28 Recently, 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.29 Here, 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 E0 state and H2 dissociation from E4 and E2.
The 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 quantum-refinement study30 and it is also supported by another study.32
The 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 133915 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.30
All MM calculations were performed with the Amber software.34 For the protein, we used the Amber ff14SB force field35 and water molecules were described by the TIP3P model.36 For the metal sites, the MM parameters were the same as in our previous investigation.30 The metal sites were treated by a non-bonded model37 and charges were obtained with the restrained electrostatic potential method, obtained at the TPSS/def2-SV(P) level of theory38,39 and sampled with the Merz–Kollman scheme.40
In 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,43 All MM atoms were included in the point-charge model, except the CL atoms.44
The total QM/MM energy in ComQum was calculated as41,42
(1) |
It should be noted that the M06-2X and M06-HF functionals were constructed for main-group thermochemistry and charge-transfer excitations, respectively, and are not recommended for systems with significant static correlation, like transition metals.56 The same applies also to BHLYP with its large amount of HF exchange.63 These 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 MoFe7S9C(homocitrate) (CH3S)(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 structure7 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 En state, i.e. the number of added protons) and it is shown in Fig. 1. Following recent Mössbauer, anomalous dispersion and QM investigations,13,24,64,65 we used the oxidation state-assignment MoIIIFeII3FeIII4 of the metal ions in the resting (E0) 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 E0, E2 and E4. One proton was added together with each electron, so the net charge of the cluster was −3 for all three states.
Fig. 1 The FeMo cluster illustrating the QM system used in all QM/MM calculations and the naming of the atoms in the cluster. |
Experiments have shown that the ground spin state of E0 and E2 are quartets with a surplus of three α electrons, whereas E4 is a doublet.3,13 These 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 states31 found in our previous study.29 These 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 collected in Table S6 (ESI†), whereas Table S7 (ESI†) list the coordinates for the complexes.
E | Position of the H atomsa | BS Stateb | # Fe–X bonds | # other bonds | n | |||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
H1 | H2 | H3 | H4 | Fe | S | C | H | S–H | C–H | H–H | ||||
a S2B(3) means that the proton on S2B points towards S3A; C(2367) means that the proton is on the Fe2–Fe3–Fe6–Fe7 face and similar for C(3457) and C(2456); S2A(Fe1) and S2A(Mo) mean that the proton on S2A points towards Fe1 or Mo, respectively; Fe2/6(5) or (3) means that the hydride ion bridging Fe2 and Fe6 points towards S5A or S3A; S5A(2) means that the proton on S5A points towards S2B; Fe3/7(2) means that the hydride ion bridging Fe3 and Fe7 points towards S2B; hydride ions terminally bound to one Fe ion are trans to the central carbide ion. b The BS states are defined by which three Fe ions have a negative spin,31viz. Fe2, 3, and 5 for BS7-1, Fe2, 3 and 4 for BS2, Fe1, 2 and 5 for BS10-1, Fe1, 3 and 7 for BS9-2, Fe1, 5 and 7 for BS6-2, Fe3, 4 and 6 for BS7-3 and Fe3, 4 and 5 for BS8-6. c H2 bound end-on to Fe5. d H2 bound side-on to Fe2. | ||||||||||||||
0 | E0 | BS7-1 | 12 | 22 | 6 | 0 | 0 | 0 | 0 | 3 | ||||
1 | E4 | S2B(3) | C(2367) | C(3457) | C(2456) | BS2 | 4 | 21 | 2 | 0 | 1 | 3 | 0 | 7 |
2 | E4 | S2B(3) | C(2367) | C(3457) | S2A(Fe1) | BS2 | 8 | 21 | 3 | 0 | 2 | 2 | 0 | 7 |
3 | E4 | S2B(3) | C(2367) | Fe6/7 | S2A(Mo) | BS10-1 | 9 | 21 | 5 | 1 | 2 | 1 | 0 | 5 |
4 | E4 | S2B(3) | Fe2/6(5) | S5A(2) | Fe3/7(2) | BS9-2 | 12 | 22 | 6 | 2 | 2 | 0 | 0 | 3 |
5 | E4 | S2B(3) | Fe2/6(3) | Fe4 | Fe5 | BS6-2 | 12 | 22 | 6 | 3 | 1 | 0 | 0 | 1 |
6 | E4 | S2B(3) | Fe2/6(3) | Fe5c | H2 | BS7-3 | 12 | 22 | 6 | 2 | 1 | 0 | 1 | 3 |
7 | E4 | S2B(3) | Fe2d | Fe4 | Fe2d | BS7-3 | 12 | 22 | 6 | 3 | 1 | 0 | 1 | 3 |
8 | E2 | C(3457) | C(2367) | BS8-6 | 7 | 22 | 4 | 0 | 0 | 2 | 0 | 5 | ||
9 | E2 | S2B(3) | Fe2/6(3) | BS7-3 | 12 | 22 | 6 | 1 | 1 | 0 | 0 | 3 |
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.
From 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.28 However, a similar study based on a GGA functional would have reached the opposite conclusion.
Fig. 3 Energies for the 1 → 5, 4 → 5, 3 → 4, 2 → 3 and 1 → 2 isomerisations, obtained with ten different DFT methods. The cyan line shows % HF for the various methods (right axis). |
Thus, the various DFT methods give ΔE15 results that differ by up to 595 kJ mol−1. This is 4–30 times larger variation than what is found in other systems.67–69 For example, Fig. 4 shows the reaction energies of the N2 + 3H2 → 2NH3 reaction, divided into three steps (involving the formation of N2H2, N2H4 and finally two NH3 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 N2H2) 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.
Fig. 4 Reaction energies for the N2 + H2 → N2H2 (top, blue), N2H2 + H2 → N2H4 (middle, red) and N2H4 + H2 → 2 NH3 (bottom, yellow) reactions, obtained in the gas phase with the 13 DFT methods. |
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 ΔE15 and % HF, R2 = 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 systems56,63 and they gave qualitative changes in the structure of 5, we present these results in the ESI.† For the N2 + 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 ΔE15 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 E4 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.27 From Fig. 3 (red line), it can be seen that ΔE45 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 ΔE45 and % HF, R2 = 0.8.
In model 3, there is only one hydride ion and the carbide ion is singly protonated. Interestingly, ΔE34 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 ΔE34 and % HF (R2 = 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 ΔE23 and % HF (R2 = 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 ΔE12 and % HF (0.3, but it increases to 0.75 with all 13 DFT functionals).
Thus, the large variation in ΔE15 among the DFT functionals is caused mainly by ΔE34 and ΔE23. 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 E0 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.
Next, we studied two structures of the E2 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 μ2 sulfide and one hydride ion bridging Fe2 and Fe6 (9; cf.Fig. 2). They are the best E2 structures obtained in our previous study with B3LYP and TPSS, respectively.29 As 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 ΔE89 to % HF is 0.6.
With these structures, we can estimate the dissociation energy of H2 from either 5 or 1. It is favourable for all GGA functionals (5 → 9 + H2) by 68–152 kJ mol−1. For the hybrid functionals (1 → 8 + H2), H2 dissociation is also favourable by 3–112 kJ mol−1, except for PBE0 (Fig. 5, green line).
Finally, we studied also the resting state of the enzyme (E0; 0) to obtain an estimate of the dissociation energy of H2 from E2 (9/8 → 0 + H2). 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.
Fig. 6 Performance of the various DFT methods for the structure of the resting E0 state compared to the crystal structure.6 For 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 (<3 Å) metal–metal (M–M) distances. |
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 (≤−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.8 on average. For structures 1 and 2, these functionals indicate a high-spin state of Mo (spin populations of 1.6–2.5), which are not consistent with 95Mo ENDOR experiments.70–72
Likewise, for the homolytic Co(III)–C bond dissociation energy in methyl corrin (coenzyme B12 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).68 Later, 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 R2 = 0.89.73
For 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.74 The 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.69
Thus, 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 main-group 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).75 Hybrid functionals remain better for some transition-metal complexes, in particular with hard ligands.67,76 However, 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 models78,79 and for complexes with nitride and carbonyl ligands.80 In 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.81 Clearly, the performance of DFT methods depend on the system studied, probably depending on the importance of static correlation.80
An 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 E4 state has two bridging hydride ions.3 Unfortunately, the results in Fig. 3 show that none of the DFT functionals suggest that the best E4 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 N2. Of course, it is then essential that H2 does not form and dissociate before N2 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 H2 dissociation energies from both E4 and E2.
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.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: 10.1039/c8cp06930a |
This journal is © the Owner Societies 2019 |