Open Access Article
This Open Access Article is licensed under a
Creative Commons Attribution 3.0 Unported Licence

Extremely large differences in DFT energies for nitrogenase models

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

Received 8th November 2018 , Accepted 20th December 2018

First published on 21st December 2018


Abstract

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.


Introduction

Nitrogenase (EC 1.18/19.6.1) is the only enzyme in nature that can cleave the triple bond in N2 to form ammonia and make nitrogen available for cell metabolism.1–3 It is present in a few groups of bacteria and archaea.1–3 The nitrogenase reaction is essential to the life on earth – although 78% of the atmosphere is N2, nitrogen is often the limiting factor for plant growth and a main component of synthetic fertilisers.3 In fact, the industrial Haber–Bosch process to form ammonia consumes ∼1% of the world's energy supplies, produces half of the total biologically available nitrogen on earth4 and is a major factor in the recent human population explosion.3

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.

Methods

The protein

All calculations are based on the 1.0 Å crystal structure of nitrogenase from Azotobacter vinelandii (PDB code 3U7Q).7 The setup of the protein is identical to that of our previous studies.29–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.30

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 133[thin space (1/6-em)]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.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

QM/MM calculations

The QM/MM calculations were performed with the ComQum software.41,42 In 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).30

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

 
image file: c8cp06930a-t1.tif(1)
where EHLQM1+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). image file: c8cp06930a-t2.tif is the MM energy of the QM system, still truncated by HL atoms, but without any electrostatic interactions. Finally, image file: c8cp06930a-t3.tif 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.45 The 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.

QM calculations

All QM calculations were performed with the Turbomole software (versions 7.1 and 7.2).46 We employed 13 different DFT methods, TPSS,38 PBE,47 BP86,48,49 BLYP,48,50 B97D,51 PBE0,52 TPSS,53 B3LYP,48,50,54 BHLYP,55 M06,56 M06-L,57 M06-2X,56 and M06-HF.58 TPSS, PBE, BP86, BLYP, B97D and M06-L are (meta-) generalized gradient approximation (GGA) functionals, with no admixture of Hartree–Fock (HF) exchange. The other DFT methods are hybrid functionals with 10% (TPSSh), 20% (B3LYP), 25% (PBE0), 27% (M06), 50% (BHLYP), 54% (M06-2X), and 100% (M06-HF) HF exchange (% HF). All calculations employed the def2-SV(P) basis set,39 because previous studies have shown that increasing the basis set to def2-TZVPD changes the relative energies by less than 16 kJ mol−1.29,31 The calculations were sped up by expanding the Coulomb interactions in an auxiliary basis set, the resolution-of-identity approximation.59,60 Empirical dispersion corrections were included with the DFT-D3 approach61 and Becke–Johnson damping,62 using parameters optimised for each DFT method, as implemented in Turbomole for all functionals except M06, M06-L, M06-2X and M06-HF (for which no parameters are available because the methods should include dispersion effects by parameterisation).

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.


image file: c8cp06930a-f1.tif
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.

Table 1 Description of the ten models, listing the En state, the positions of the added protons (atom names are shown in Fig. 1; the number in brackets indicates the direction of the protona[thin space (1/6-em)]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)
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.

Results

In this study, we compare energies and geometries of 10 models of the FeMo cluster in the E4, E2 or E0 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.
image file: c8cp06930a-f2.tif
Fig. 2 Structures of the 10 studied nitrogenase models, 09. 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.

E4 states

We first compare the relative stability of the best protonation states of E4 obtained with either the B3LYP (1) or TPSS (5) functionals in our extensive study, ΔE15.29 As can be seen from Fig. 2, both structures are protonated on one of the μ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.66

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.


image file: c8cp06930a-f3.tif
Fig. 3 Energies for the 15, 45, 34, 23 and 12 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.


image file: c8cp06930a-f4.tif
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 45, 34 and 23 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, 34, 23 and 12 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.

H2 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 H2 from the E4 and E2 states of the MoFe cluster. First, we included two complexes with H2 bound to the FeMo cluster in the E4 state, both with a proton on S2B. Structure 6 has H2 bound end-on to Fe5 and the fourth H is bridging Fe2 and Fe6, whereas in structure 7, H2 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 E4 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.
image file: c8cp06930a-f5.tif
Fig. 5 Energies for the 67, 5/16, 89, 5/19/8 + H2 and 9/80 + H2 reactions, obtained with ten different DFT methods (using the structure with the lowest energy for 5/1 and 9/8 for each method).

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 (59 + H2) by 68–152 kJ mol−1. For the hybrid functionals (18 + 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/80 + 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.

Geometry of the resting state

For E0, we can also compare the QM/MM structures with the crystal structure of nitrogenase (3U7Q).7 It 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 broken-symmetry 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.
image file: c8cp06930a-f6.tif
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.

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).29 These 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 time-consuming 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 (≤−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

Discussion

We have studied how the relative energies of various protonation and H2-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)S4NL complexes with L = NH3, N2H4, CO, NO+, PH3 and PMe3.67 They 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.28

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.

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,28 Unfortunately, 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 H2-dissociation energies but a high-spin Mo state for the best E4 structures. None of the tested methods indicate that the best E4 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.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work has been supported by grants from the Swedish Research Council (project 2014-5540), from COST through Action CM1305 (ECOSTBio), by the Royal Physiographic Society in Lund and by a scholarship to LC from the China Scholarship Council. The computations were performed on computer resources provided by the Swedish National Infrastructure for Computing (SNIC) at Lunarc at Lund University and HPC2N at Umeå University.

References

  1. B. K. Burgess and D. J. Lowe, Chem. Rev., 1996, 96, 2983–3012 CrossRef CAS PubMed.
  2. B. Schmid, H.-J. Chiu, V. Ramakrishnan, J. B. Howard and D. C. Rees, in Handbook of Metalloproteins, John Wiley & Sons, Ltd, 2006, pp. 1025–1036 Search PubMed.
  3. B. M. Hoffman, D. Lukoyanov, Z.-Y. Yang, D. R. Dean and L. C. Seefeldt, Chem. Rev., 2014, 114, 4041–4062 CrossRef CAS PubMed.
  4. B. E. Smith, Science, 2002, 297, 1654–1655 CrossRef CAS PubMed.
  5. J. Kim and D. C. Rees, Science, 1992, 257, 1677–1682 CrossRef CAS PubMed.
  6. O. Einsle, F. A. Tezcan, S. L. A. Andrade, B. Schmid, M. Yoshida, J. B. Howard and D. C. Rees, Science, 2002, 297, 1696 CrossRef CAS PubMed.
  7. T. Spatzal, M. Aksoyoglu, L. Zhang, S. L. A. Andrade, E. Schleicher, S. Weber, D. C. Rees and O. Einsle, Science, 2011, 334, 940 CrossRef CAS PubMed.
  8. T. Spatzal, K. A. Perez, O. Einsle, J. B. Howard and D. C. Rees, Science, 2014, 345, 1620–1623 CrossRef CAS PubMed.
  9. O. Einsle, J. Biol. Inorg. Chem., 2014, 19, 737–745 CrossRef CAS PubMed.
  10. R. R. Eady, Chem. Rev., 1996, 96, 3013–3030 CrossRef CAS PubMed.
  11. R. N. F. Thorneley and D. J. Lowe, in Molybdenum Enzymes, ed. T. G. Spiro, Wiley, New York, 1985, pp. 221–284 Search PubMed.
  12. K. M. Lancaster, M. Roemelt, P. Ettenhuber, Y. Hu, M. W. Ribbe, F. Neese, U. Bergmann and S. DeBeer, Science, 2011, 334, 974–977 CrossRef CAS PubMed.
  13. R. Bjornsson, F. A. Lima, T. Spatzal, T. Weyhermüller, P. Glatzel, E. Bill, O. Einsle, F. Neese and S. DeBeer, Chem. Sci., 2014, 5, 3096–3103 RSC.
  14. F. Tuczek, in RSC Metallobiology Series 7, ed. R. Hille, C. Schulzke and M. L. Kirk, Royal Society of Chemistry, Cambridge, 2017, pp. 223–274 Search PubMed.
  15. I. Dance, J. Biol. Inorg. Chem., 1996, 1, 581–586 CrossRef CAS.
  16. K. K. Stavrev and M. C. Zerner, Int. J. Quantum Chem., 1998, 70, 1159–1168 CrossRef CAS.
  17. P. E. M. Siegbahn, J. Westerberg, M. Svensson and R. H. Crabtree, J. Phys. Chem. B, 1998, 102, 1615–1623 CrossRef CAS.
  18. T. Lovell, J. Li, T. Liu, D. A. Case and L. Noodleman, J. Am. Chem. Soc., 2001, 123, 12392–12410 CrossRef CAS PubMed.
  19. H. Xie, R. Wu, Z. Zhou and Z. Cao, J. Phys. Chem. B, 2008, 112, 11435–11439 CrossRef CAS PubMed.
  20. J. Kästner and P. E. Blöchl, J. Am. Chem. Soc., 2007, 129, 2998–3006 CrossRef PubMed.
  21. P. P. Hallmen and J. Kästner, Z. Anorg. Allg. Chem., 2015, 641, 118–122 CrossRef CAS.
  22. I. Dance, Z. Anorg. Allg. Chem., 2015, 641, 91–99 CrossRef CAS.
  23. J. B. Varley, Y. Wang, K. Chan, F. Studt and J. K. Nørskov, Phys. Chem. Chem. Phys., 2015, 17, 29541–29547 RSC.
  24. P. E. M. Siegbahn, J. Am. Chem. Soc., 2016, 138, 10485–10495 CrossRef CAS PubMed.
  25. M. L. McKee, J. Phys. Chem. A, 2016, 120, 754–764 CrossRef CAS PubMed.
  26. L. Rao, X. Xu and C. Adamo, ACS Catal., 2016, 6, 1567–1577 CrossRef CAS.
  27. D. Lukoyanov, N. Khadka, D. R. Dean, S. Raugei, L. C. Seefeldt and B. M. Hoffman, Inorg. Chem., 2017, 56, 2233–2240 CrossRef CAS PubMed.
  28. P. E. M. Siegbahn, J. Comput. Chem., 2018, 39, 743–747 CrossRef CAS PubMed.
  29. L. Cao and U. Ryde, J. Chem. Theory Comput., 2018, 14, 6653–6678 CrossRef CAS PubMed.
  30. L. Cao, O. Caldararu and U. Ryde, J. Phys. Chem. B, 2017, 121, 8242–8262 CrossRef CAS PubMed.
  31. L. Cao and U. Ryde, Int. J. Quantum Chem., 2018, 118, e25627 CrossRef.
  32. B. Benediktsson and R. Bjornsson, Inorg. Chem., 2017, 56, 13417–13429 CrossRef CAS PubMed.
  33. B. M. Barney, J. Mcclead, D. Lukoyanov, M. Laryukhin, T. Yang, D. R. Dean, B. M. Hoffman and L. C. Seefeldt, Biochemistry, 2007, 46, 6784–6794 CrossRef CAS PubMed.
  34. D. A. Case, J. T. Berryman, R. M. Betz, D. S. Cerutti, T. E. Cheatham, T. A. Darden, R. E. Duke, T. J. Giese, H. Gohlke, A. W. Goetz, N. Homeyer, S. Izadi, P. Janowski, J. Kaus, A. Kovalenko, T. S. Lee, S. LeGrand, P. Li, T. Luchko, R. Luo, B. Madej, K. M. Merz, G. Monard, P. Needham, H. Nguyen, H. T. Nguyen, I. Omelyan, A. Onufriev, D. R. Roe, A. E. Roitberg, R. Salomon-Ferrer, C. Simmerling, W. Smith, J. Swails, R. C. Walker, J. Wang, R. M. Wolf, X. Wu, D. M. York and P. A. Kollman, AMBER 2014, University of California, San Francisco, 2014.
  35. J. A. Maier, C. Martinez, K. Kasavajhala, L. Wickstrom, K. E. Hauser and C. Simmerling, J. Chem. Theory Comput., 2015, 11, 3696–3713 CrossRef CAS PubMed.
  36. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, J. Chem. Phys., 1983, 79, 926–935 CrossRef CAS.
  37. L. Hu and U. Ryde, J. Chem. Theory Comput., 2011, 7, 2452–2463 CrossRef CAS PubMed.
  38. J. Tao, J. P. Perdew, V. N. Staroverov and G. E. Scuseria, Phys. Rev. Lett., 2003, 91, 146401 CrossRef PubMed.
  39. A. Schäfer, H. Horn and R. Ahlrichs, J. Chem. Phys., 1992, 97, 2571–2577 CrossRef.
  40. B. H. Besler, K. M. Merz and P. A. Kollman, J. Comput. Chem., 1990, 11, 431–439 CrossRef CAS.
  41. U. Ryde, J. Comput.-Aided Mol. Des., 1996, 10, 153–164 CrossRef CAS PubMed.
  42. U. Ryde and M. H. M. Olsson, Int. J. Quantum Chem., 2001, 81, 335–347 CrossRef CAS.
  43. N. Reuter, A. Dejaegere, B. Maigret and M. Karplus, J. Phys. Chem. A, 2000, 104, 1720–1735 CrossRef CAS.
  44. L. Hu, P. Söderhjelm and U. Ryde, J. Chem. Theory Comput., 2011, 7, 761–777 CrossRef CAS PubMed.
  45. L. Cao and U. Ryde, Front. Chem., 2018, 6, 89 CrossRef PubMed.
  46. F. Furche, R. Ahlrichs, C. Hättig, W. Klopper, M. Sierka and F. Weigend, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2014, 4, 91–100 CAS.
  47. J. P. Perdew, K. Burke, M. Ernzerhof, B. Kieron and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868 CrossRef CAS PubMed.
  48. A. D. Becke, Phys. Rev. A, 1988, 38, 3098–3100 CrossRef CAS.
  49. J. P. Perdew, Phys. Rev. B, 1986, 33, 8822–8824 CrossRef.
  50. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B, 1988, 37, 785–789 CrossRef CAS.
  51. S. Grimme, J. Comput. Chem., 2006, 27, 1787–1799 CrossRef CAS PubMed.
  52. J. P. Perdew, M. Ernzerhof and K. Burke, J. Chem. Phys., 1996, 105, 9982 CrossRef CAS.
  53. V. N. Staroverov, G. E. Scuseria, J. Tao and J. P. Perdew, J. Chem. Phys., 2003, 119, 12129–12137 CrossRef CAS.
  54. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652 CrossRef CAS.
  55. A. D. Becke, J. Chem. Phys., 1993, 98, 1372 CrossRef CAS.
  56. Y. Zhao and D. G. Truhlar, Theor. Chem. Acc., 2008, 120, 215–241 Search PubMed.
  57. Y. Zhao and D. G. Truhlar, J. Chem. Phys., 2006, 125, 194101 CrossRef PubMed.
  58. Y. Zhao and D. G. Truhlar, J. Phys. Chem. A, 2006, 110, 13126–13130 CrossRef CAS PubMed.
  59. K. Eichkorn, O. Treutler, H. Öhm, M. Häser and R. Ahlrichs, Chem. Phys. Lett., 1995, 240, 283–289 CrossRef CAS.
  60. K. Eichkorn, F. Weigend, O. Treutler and R. Ahlrichs, Theor. Chem. Acc., 1997, 97, 119–124 Search PubMed.
  61. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104 CrossRef PubMed.
  62. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465 CrossRef CAS PubMed.
  63. H. S. Yu, S. L. Li and D. G. Truhlar, J. Chem. Phys., 2016, 145, 130901 CrossRef PubMed.
  64. T. Spatzal, J. Schlesier, E.-M. Burger, D. Sippel, L. Zhang, S. L. A. Andrade, D. C. Rees and O. Einsle, Nat. Commun., 2016, 7, 10902 CrossRef CAS PubMed.
  65. R. Bjornsson, F. Neese and S. DeBeer, Inorg. Chem., 2017, 56, 1470–1477 CrossRef CAS PubMed.
  66. P. E. M. Siegbahn, Inorg. Chem., 2018, 57, 1090–1095 CrossRef CAS PubMed.
  67. M. Reiher, O. Salomon and B. A. Hess, Theor. Chem. Acc., 2001, 107, 48–55 Search PubMed.
  68. K. P. Jensen and U. Ryde, J. Phys. Chem. A, 2003, 107, 7539–7545 CrossRef CAS.
  69. J. Li, M. Andrejic, R. A. Mata and U. Ryde, Eur. J. Inorg. Chem., 2015, 3580–3589 CrossRef CAS.
  70. S. Raugei, L. C. Seefeldt and B. M. Hoffman, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 10521–10530 CrossRef PubMed.
  71. R. A. Venters, M. J. Nelson, P. A. McLean, A. E. True, M. A. Levy, B. M. Hoffman and W. H. Orme-Johnson, J. Am. Chem. Soc., 1986, 108, 3487–3498 CrossRef CAS.
  72. D. Lukoyanov, Z.-Y. Yang, D. R. Dean, L. C. Seefeldt and B. M. Hoffman, J. Am. Chem. Soc., 2010, 132, 2526–2527 CrossRef CAS PubMed.
  73. U. Ryde, R. A. Mata and S. Grimme, Dalton Trans., 2011, 40, 11176–11183 RSC.
  74. J. Li, R. A. Mata and U. Ryde, J. Chem. Theory Comput., 2013, 9, 1799–1807 CrossRef CAS PubMed.
  75. L. Goerigk, A. Hansen, C. Bauer, S. Ehrlich, A. Najibi and S. Grimme, Phys. Chem. Chem. Phys., 2017, 19, 32184–32215 RSC.
  76. M. Lundberg and P. E. M. Siegbahn, J. Comput. Chem., 2005, 26, 661–667 CrossRef CAS PubMed.
  77. I. Dance, Mol. Simul., 2018, 44, 568–581 CrossRef CAS.
  78. M. G. Delcey, K. Pierloot, Q. M. Phung, S. Vancoillie, R. Lindh and U. Ryde, Phys. Chem. Chem. Phys., 2014, 16, 7927–7938 RSC.
  79. G. Dong, Q. M. Phung, S. D. Hallaert, K. Pierloot and U. Ryde, Phys. Chem. Chem. Phys., 2017, 19, 10590–10601 RSC.
  80. J. N. Harvey, Annu. Rep. Prog. Chem., Sect. C: Phys. Chem., 2006, 102, 203–226 RSC.
  81. S. Dohm, A. Hansen, M. Steinmetz, S. Grimme and M. P. Checinski, J. Chem. Theory Comput., 2018, 14, 2596–2608 CrossRef CAS PubMed.

Footnote

Electronic supplementary information (ESI) available. See DOI: 10.1039/c8cp06930a

This journal is © the Owner Societies 2019
Click here to see how this site uses Cookies. View our privacy policy here.