What Is the Structure of the E4 Intermediate in Nitrogenase?

Nitrogenase is the only enzyme that can cleave the strong triple bond in N2. The active site contains a complicated MoFe7S9C cluster. It is believed that it needs to accept four protons and electrons, forming the E4 state, before it can bind N2. However, there is no consensus on the atomic structure of the E4 state. Experimental studies indicate that it should contain two hydride ions bridging two pairs of Fe ions, and it has been suggested that both hydride ions as well as the two protons bind on the same face of the cluster. On the other hand, density functional theory (DFT) studies have indicated that it is energetically more favorable with either three hydride ions or with a triply protonated carbide ion, depending on the DFT functional. We have performed a systematic combined quantum mechanical and molecular mechanical (QM/MM) study of possible E4 states with two bridging hydride ions. Our calculations suggest that the most favorable structure has hydride ions bridging the Fe2/6 and Fe3/7 ion pairs. In fact, such structures are 14 kJ/mol more stable than structures with three hydride ions, showing that pure DFT functionals give energetically most favorable structures in agreement with experiments. An important reason for this finding is that we have identified a new type of broken-symmetry state that involves only two Fe ions with minority spin, in contrast to the previously studied states with three Fe ions with minority spin. The energetically best structures have the two hydride ions on different faces of the FeMo cluster, whereas better agreement with ENDOR data is obtained if they are on the same face; such structures are only 6–22 kJ/mol less stable.


■ INTRODUCTION
Nitrogenase (EC 1.18/19.6.1) is the only enzyme that can cleave the triple bond of N 2 , forming two molecules of ammonia and making atmospheric nitrogen available for biological systems. 1−3 Nitrogenase is found only in some bacteria and cyanobacteria, but several higher plants, e.g. legumes, rice, and alder, live in symbiosis with such microorganisms. Nitrogenase consists of two proteins, the Fe protein that delivers electrons and the MoFe protein that contains the active site. Crystallographic studies have shown that the Fe protein contains a Fe 4 S 4 cluster, whereas the MoFe protein contains the Fe 8 S 7 Cys 6 Pcluster, used for electron transfer, and the catalytic MoFe 7 S 9 C-(homocitrate)CysHis FeMo cluster (shown in Figure 1). 4−8 Alternative nitrogenases exist, in which the Mo ion is replaced by V or Fe. 9 Nitrogenases catalyze the chemical reaction The binding of two molecules of ATP triggers the binding of the reduced Fe protein to the MoFe protein and the delivery of one electron to the FeMo cluster. 1−3 Then, both ATP molecules are hydrolyzed, and the Fe protein dissociates. This is repeated eight times during the catalytic cycle. Consequently, the reaction mechanism of nitrogenase is normally described by nine intermediates, E 0 −E 8 , differing in the number of delivered electrons, the Lowe−Thorneley cycle. 10 E 0 is the resting state, which is in the Mo III Fe 3 II Fe 4 III oxidation state. 11−14 It is currently believed that four electrons and protons need to be added to E 0 , i.e., forming the E 4 state, before N 2 may bind. 1−3,12,15 N 2 binding is facilitated by the reductive elimination of H 2 , which would explain why H 2 is a necessary byproduct of the reaction. 3 Many computational studies of nitrogenase have been presented, with the hope to give a detailed atomistic description of the various reaction steps. 3,12,23−28,13,16−22 Unfortunately, they have given strongly divergent results. For example, some studies have suggested that N 2 is protonated alternatively on both N atoms, giving diazene (HNNH) and hydrazine (H 2 NNH 2 ) as intermediates, 24,27 whereas others have suggested that the distal N atom is first triply protonated and dissociates as NH 3 before the second N atom is protonated. 25,26 Likewise, some studies have suggested that N 2 binds end-on to one 24 or two Fe ions, 25,28 whereas others have suggested that it binds side-on to two Fe ions 28 or that it forms a covalent bond with the central carbide ion. 26,27 In fact, there is not any consensus regarding the structure of the central E 4 intermediate that binds N 2 . It should contain four protons, and Hoffman and co-workers have suggested, based on ENDOR experiments, that two of these bind to sulfide ions, whereas the other two bridge between two pairs of Fe ions as hydride ions. 3,29,30 With the help of density functional theory (DFT) calculations, they have suggested a structure in which both hydride ions and the two protons bind on the same face of the cluster, viz., between the Fe2/6 and Fe3/7 pairs, as well as on the belt S2B and S5A ions, as is shown in Figure 2 (the name of the Fe and S ions, taken from the crystal structure, 6 are also shown in the figure). 28,30 On the other hand, Siegbahn has argued that it is energetically much more favorable if the protons bind to the central carbide, so that the most stable structure with four electrons and protons has one proton on a sulfide ion and three on the carbide, forming CH 3 which moves out from the center of the cluster (giving rise to a cavity where N 2 can subsequently bind). 13,31,32 We have recently shown that the reason for this discrepancy is that different DFT methods give very different results for such protonation isomers. 33,34 Pure functionals, like BP86, used by Hoffman and co-workers, strongly favor hydrogen atoms binding to Fe ions as hydride ions, whereas hybrid functionals (like B3LYP, used by Siegbahn), favor proton binding to sulfide ions and especially to the central carbide ion. In fact, predictions of the relative energies of the various protonation states may differ by up to 600 kJ/mol, depending on the functional. Moreover, finding the most stable protonation state of E 4 is a formidable task. There are at least 50 different possible positions and conformations of each proton, 33,35 which together with the uncertainty in the spin and broken-symmetry states give over 700 million possibilities for the E 4 state. 33 Using a heuristic, but systematic approach, we have suggested optimum protonation states for the E 0 −E 4 intermediates, obtained with the TPSS and B3LYP methods. 33 For E 4 , TPSS suggested structures with three hydride ions and one proton on a sulfide ion. Many structures were close in energy, and the hydride ions seem to be able to move rather freely between different Fe ions. On the other hand, B3LYP strongly disfavors metal-bound hydride ions and instead suggests that E 4 involves a triply protonated carbide and a protonated sulfide ion, in agreement with the suggestion of Siegbahn. 13,31 Unfortunately, neither of these structures is in accordance with the ENDOR experiments, suggesting instead that E 4 involves two bridging hydride ions. 3,29 A more detailed study with 13 DFT functionals confirmed this observation. None of the tested functionals pointed out Hoffman's structure with two bridging hydride ions as the most stable E 4 structure. 34 A possible explanation to this problem is that the Hoffman structure is not the best E 4 structure with two bridging hydridesit was based on only a restricted search on one of the three faces of the FeMo cluster. In this paper, we therefore perform a more systematic search for the best E 4 structure with two bridging hydride ions and two protonated sulfides, using combined QM/MM calculations 36 and methods to deal with the broken-symmetry states and other computational difficulties developed in our previous studies. 33,34,37 ■ METHODS The Protein. All calculations were based on the 1.0-Å crystal structure of nitrogenase from Azotobacter vinelandii (PDB code 3U7Q). 6 The setup of the protein is identical to that of our previous studies of the protein. 33,34,37,38 The entire heterotetramer was included in the calculations because the various subunits are entangled without any natural way to separate them. The QM calculations were concentrated on the FeMo clusters in the C subunit because there is a buried imidazole molecule from the solvent rather close to the active site (∼11 Å) in the A subunit. The P clusters and the FeMo cluster in subunit A were modeled by MM in the fully reduced and resting states, respectively. 38 The protonation states of all residues were the same as before. 38 All Arg, Lys, Asp, and Glu residues were assumed to be 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 residues). Cys residues coordinating to Fe ions were assumed to be 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 presumed to be protonated on both the ND1 and NE2 atoms (and therefore positively charged), whereas the remaining 14 His residues were modeled with a proton on the NE2 atom. The Structure of E 4 suggested by Hoffman and co-workers, 3,28 as obtained by our QM/MM geometry optimizations. The protons (shown as green balls) bind to the S2B and S5A ions, whereas the hydride ions (also green balls) bridge the Fe2/6 and Fe3/7 pairs. It represents the S2B(5)−S5A(2)−Fe2/6(5)−Fe3/7(2) conformation in our nomenclature. The figure also shows the names of the Fe and S ions, following the 3U7Q crystal structure. 6 homocitrate was modeled 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 one in a recent extensive QM/MM, molecular dynamics, and quantum-refinement study, 38 and this protonation state is also supported by another recent study. 14 The protein was solvated in a sphere with a radius of 65 Å around the geometrical center of the protein. 160 Cl − and 182 Na + ions were added at random positions to neutralize the protein and give an ionic strength of 0.2 M. 39 The final system contained 133 915 atoms. The added protons, counterions, and water molecules were optimized by a simulated annealing calculation (up to 370 K), followed by a minimization, keeping the other atoms fixed at the crystal-structure positions. 38 All MM calculations were performed with the Amber software. 40 For the protein, we used the Amber ff14SB force field, 41 and water molecules were described by the TIP3P model. 42 For the metal sites, the MM parameters were the same as in our previous investigation. 38 The only exception is the Ca ion at the subunit interface, which was replaced by a Fe ion, in accordance with recent crystallographic studies. 43 The metal sites were treated by a nonbonded model, 44 and charges were obtained with the restrained electrostatic potential method, obtained at the TPSS/def2-SV(P) level of theory 45,46 and sampled with the Merz−Kollman scheme. 47 The charges of the Fe site are listed in Table S1 in the Supporting Information QM Calculations. All QM calculations were performed with the Turbomole software (versions 7.1 and 7.2). 48 Most calculations employed the TPSS 45 method and the def2-SV(P) 46 basis set. However, in some calculations, we used the much larger def2-TZVPD 49 basis set, although previous studies have shown that the energies typically change only by up to 11− 15 kJ/mol when the basis set is enlarged. 33,34,37 Moreover, in some calculations, we employed seven additional DFT methods: PBE, 50 M06-L, 51 B97D, 52 PBE0, 53 TPSSh, 54 B3LYP, 55−57 and M06. 58 All DFT calculations were sped up by expanding the Coulomb interactions in an auxiliary basis set, the resolution-ofidentity (RI) approximation. 59,60 Empirical dispersion corrections were included with the DFT-D3 approach 61 and Becke− Johnson damping, 62 as implemented in Turbomole.
The FeMo cluster was modeled by MoFe 7 S 9 C(homocitrate)-(CH 3 S)(imidazole), where the two last groups are models of Cys-275 and His-442. In addition, all groups that form hydrogen bonds to the FeMo cluster in the crystal structure 6 were also included, viz., Arg-96 and His-195 (side chains); Ser-278 and Arg-359 (both backbone and side chain, including the Cα and C and O atoms from Arg-277); and Gly-356, Gly-357, and Leu-358 (backbone, including the Cα and C and O atoms from Ile-355), as well as two water molecules, in total 151 atoms (shown in Figure 1). For the best structures, we also tested a larger QM system, including also the side chains of Val-70 and Gln-191, the full side chain of Arg-96, and an extra methyl group for the two His residues (186 atoms in total, shown in Figure S1). Following recent Mossbauer, anomalous dispersion, and QM investigations, 11 ). 37 The various BS states are named by listing the number in the Noodleman nomenclature (BS1-10), 20 followed by the numbers of the three Fe ions with minority spin, e.g., BS7-235, indicating that Fe2, Fe3, and Fe5 have β spin. This is slightly different from what was used in our previous studies 33,34,37 (but more transparent) and is the same used by Benediktsson and Bjornsson. 63 We provide a translation between the two nomenclatures in Table S2 in the Supporting  Information. A starting wave function for the FeMo cluster was obtained by first optimizing the all-high-spin state with 35 unpaired electrons and then changing the total α and β occupation numbers to the desired net spin. This gave one of the BS states. The other BS states were then obtained by simply swapping the coordinates of the Fe ions. 64 In many cases, we instead used the fragment approach by Szilagyi and Winslow to obtain a proper BS state. 65 As is discussed below, binding of hydride ions reduces the spin population on the corresponding Fe ion and makes the electronic states more complicated. In fact, we found in several cases that a new type of BS state with only two Fe ions of β spin (and therefore five with α spin) was most favorable. Such states can be obtained in 21 different ways ( 7 6 2 × ), and they are denoted BS-ij, where i and j are the two Fe ions with minority spin, e.g., .
We have previously studied the 35 BS states for the resting state, the protonated resting state, and the reduced state, and how their energies vary with the QM method, the size of the basis set, the geometry, and the influence of the surroundings. 37 The conclusions were that the effects of the basis set and the surroundings were restricted (up to 7−11 kJ/mol), and the effect of the geometry was intermediate (up to 37 kJ/mol, but the correlation, R 2 , was 0.92−0.98). Also, the effect of the DFT functional (TPSS or B3LYP) was large (up to 58 kJ/mol). Therefore, we first studied all systems with the same BS state. However, in some cases, the BS state changed during the geometry optimization (some of the spin densities are often quite low and easily change sign). For the best protonation states, we performed a complete investigation of the relative energy of all 35 or 21 BS states with full geometry optimization for each state.
QM/MM Calculations. The QM/MM calculations were performed with the ComQum software. 66,67 In this approach, the protein and solvent are split into three subsystems. System 1 (the QM region) was relaxed by QM methods. System 2 contained all residues or water molecules with at least one atom within 6 Å of an atom in the QM region, in total 1377 atoms (87 residues and 35 water molecules). System 3 contained the remaining part of the protein and the solvent, and it was kept fixed at the original coordinates (equilibrated crystal structure). The total system was spherical and nonperiodic with 133 915 atoms. 38 Two sets of calculations were performed. In one, all atoms outside the QM region were kept fixed at their starting (crystallographic) positions. In the second set, atoms in system 2 were relaxed by a full MM geometry optimization (keeping systems 1 and 3 fixed) in each iteration of the QM/MM geometry optimization. The first set ensures that all energies are comparable and not affected by uncertainties in the MM force field but may present a too rigid surrounding protein. The second set samples the effects of a relaxed protein, but the structures may be affected by the risk that different structures end up in different local minima of the surrounding protein.
In the QM calculations, system 1 was represented by a wave function, whereas all the other atoms were represented by an array of partial point charges, one for each atom, taken from the MM setup. Thereby, the polarization of the QM system by the surroundings is included in a self-consistent manner (electrostatic embedding). When there is a bond between systems 1 and 2 (a junction), the hydrogen link-atom approach was employed. The QM system was capped with hydrogen atoms (hydrogen link atoms, HL), the positions of which are linearly related to the corresponding carbon atoms (carbon link atoms, CL) in the full system. 66,68 All atoms were included in the point-charge model, except the CL atoms. 69 The total QM/MM energy in ComQum was calculated as 66 where E QM1+ptch23 HL is the QM energy of the QM system truncated by HL atoms and embedded in the set of point charges modeling systems 2 and 3 (but excluding the self-energy of the point charges). E MM1,q 1 = 0 HL is the MM energy of the QM system, still truncated by HL atoms, but without any electrostatic interactions. Finally, E MM123,q 1 = 0 CL 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. 70 No cutoff was used for any of the interactions in the three energy terms in eq 3. The geometry optimizations were continued until the energy change between two iterations was less than 2.6 J/mol (10 −6 a.u.), and the maximum norm of the Cartesian gradients was below 10 −3 a.u.

■ RESULT AND DISCUSSION
In this paper, we investigate whether it is possible to find a structure of the E 4 state in nitrogenase that is consistent with the experimental ENDOR data, 3,29 i.e., with two bridging hydride ions, and that is lower in energy than other possible structures. As discussed in the Introduction and our previous study, 33 this is a formidable task because there are on the order of 10 9 possible protonation and electronic states of E 4 . Consequently, we have to make some assumptions and simplifications to make the study feasible. Thus, we first decide the best positions of the two protons on sulfide ions, assuming that the two hydride ions are at the positions suggested by Hoffman and co-workers (shown in Figure 2). 3,28 Next, we fix the positions of the two protons and test all possible structures with two bridging hydride ions. These two studies are performed with the TPSS-D3 DFT method. Finally, for the best structures, we make some variations in the positions of the protons and also in the methods. We also check the results employing seven additional DFT functionals. The results are described in separate sections below.
Protonation of Sulfide Ions. Our first aim is to find the optimum positions of the two protons in a E 4 structure with two hydride ions. In our previous study, both TPSS and B3LYP calculations agreed that the first proton (for E 1 ) binds to S2B (the names of the various atoms in the FeMo cluster, taken from the 3U7Q crystal structure, 6 are shown in Figure 2), 33 and several other studies have suggested protonation of the same atom. 13,25,35 Therefore, we always kept this atom protonated. For the second proton, we compared the energy of all possible structures obtained by protonating a sulfide ion, keeping the two bridging hydride ions on Fe2/6 and Fe3/7, as suggested by Hoffman et al. (Figure 2). 3,28 The result of this investigation is shown in Table 1. We tested 16 different structures, i.e., two different conformations for each of the remaining eight sulfide ions. For the two μ 2 (belt) sulfides, bridging the two subclusters of the FeMo cluster (S3A and S5A), the proton can point toward one of the two other belt sulfides, and this is indicated by the number of the latter belt sulfide (for example, S3A (2) indicates that the proton is pointing toward S2B, whereas S3A(5) indicates that it points toward S5A). For the μ 3 sulfides within the two subclusters, the number in parentheses indicates toward which Fe ion (or other atom, homocitrate, HCA, or water) the proton is pointing. In two The first proton binds to S2B(5), and the third and fourth hydrogen atoms bind as bridging hydrides to Fe2/6(5) and Fe3/7(2), as shown in Figure 2. The table shows the S−H or Fe−H distances (Å) of the four hydrogen atoms, as well as the relative energy (ΔE in kJ/mol), calculated either at the TPSS-D3/def2-SV(P), TPSS-D3/def2-TZVPD, or B3LYP-D3/def2-SV(P) levels of theory (the latter two single-point energies on the TPSS/def2-SV(P) structures). All complexes were in the BS4−356 state, unless otherwise stated. b This is the distance to the sulfide ion shown in the first column. c Fe7 instead of Fe2. d A terminal hydride ion, not directed toward the other Fe ion, but instead binding trans to the central carbide ion. e Changed to the BS10−147 state. f Changed to the BS3−124 state.
Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article cases, both calculations converged to the same structure, so that only 14 structures are presented in Table 1.
In one case, S3B(S1B), the binding of one of the hydride ions changed from Fe2/6 to Fe6/7. This also gave rise to a slightly changed S2B−H distance, 1.38 Å, whereas this distance was 1.36 Å for all the other structures. This structure is also the only one for which the two Fe−H distances for the bridging hydride ion are equal (1. 67  The relative energies of the various structures are shown in the last columns of Table 1. They were evaluated at the TPSS-D3/ def2-SV(P) level of theory for all structures. It can be seen that protonation of the belt S5A ion in the direction of S3A is most favorable. This structure is shown in Figure 3. Protonation of the same ion, but in the opposite direction (toward S2B), is 7 kJ/ mol less stable. The latter is the structure suggested by Hoffman and co-workers, 3,28 shown in Figure 2. At first, the stability of the S5A(2) conformation might be unexpected, as it would interfere with the hydrogen bond between the guanidinium group of Arg-96 and S5A, observed in crystal structures of the resting state (shown in Figure 1). However, in all structures in Table 1, this hydrogen bond is replaced by a strong dihydrogen bond between the same H atom of Arg-96 and the hydride ion bridging Fe3 and Fe7, with a H−H distance of only approximately 1.4 Å, shown in Figure S2.
Protonation of S3B in the direction of S1A (which led to the change from Fe2/6 to Fe6/7, discussed above) has the same energy, but it does not have two bridging hydrides. The hydride bound to Fe7 is terminal and not directed toward Fe3. Instead, it is trans to the central carbide ion, as was observed for the most stable terminal hydrides in our previous study. 33 The other structures are appreciably higher in energy (37 kJ/mol for S1A(Fe1) and 54−121 kJ/mol for the other structures).
For the three best structures, we calculated the relative energies also with the larger def2-TZVPD basis set. This changed the relative energies by only 2−3 kJ/mol, and the S5A(3) structure remained most stable. We also calculated the energies by the B3LYP-D3/def2-SV(P) (single-point energies). This changed the relative energies of the two S5A structures by only 1 kJ/mol (S5A(3) is still most stable). However, the S3B(S1A) structure was destabilized by 11 kJ/mol, as could be expected because B3LYP disfavors structures with Fe-bound hydride ions.
Thus, we conclude that it is most favorable if the second proton binds to the S5A ion in the direction toward S3A, and this structure was used in the calculations in the next section.
Structures with Two Bridging Hydride Ions. Next, we try to find the best E 4 structure with two bridging hydride ions, keeping the remaining two protons on S2B(5) and S5A (3). 95 Mo ENDOR experiments have shown that Mo does not bind the hydride ions. 71 Moreover, our previous study of the protonation states of E 0 −E 4 showed that bridging structures involving the terminal Mo or Fe1 ions were high in energy. 33 Therefore, we concentrate on structures involving the other six Fe ions. We then have three pairs of Fe ions from each subcluster (Fe2/3, Fe2/4, and Fe3/4, as well as Fe5/6, Fe5/7, and Fe6/7) and three pairs bridging the two subclusters (Fe2/6, Fe3/7, and Fe4/5). For the former six, only one conformation is possible, but for the latter, two conformations are possible, depending on the side of the bridging sulfide ion. As in our previous study, 33 we give them names indicating toward which sulfide ion they are directed. For example, Fe2/6(3) has the hydride ion on the side of S2B that is directed toward S3A, whereas Fe2/6(5) is on the opposite side, directed toward S5A. Thus, there are 12 possible positions of one bridging hydride ion, giving a total of 66 possible structures for two bridging hydride ions ( 12 11 2 × ), which were all explicitly built and optimized by QM/MM.
The results are collected in Table 2. In four cases, the hydride coordination changed to another mode during the geometry optimizations (Fe3/7(2)−Fe2/4, Fe2/3−Fe2/4, Fe2/4−Fe3/ 4, and Fe2/4−Fe5/6; these are excluded from Table 2). In all four cases, the Fe2/4 bridge changed to Fe2/6(3), which is among the most favorable mode, both in the present and in our previous study. 33 In four other structures, the two hydride ions came so close that they reacted and formed H 2 with a H−H bond length of 0.74−0.96 Å (it is 0.76 Å in free H 2 calculated at the same level of theory). This often gave rise to strongly distorted structures. In seven cases, the H−H distance was longer, 1.07−1.37 Å, giving rise to less distorted structures.
The distance between the Fe ions and the hydride ions varied between 1.47 and 2.52 Å. Sometimes, the two Fe−H bonds of the same hydride ion were almost equal (e.g., 1.58 and 1.61 Å for Fe2/6(3)−Fe2/3), but it was more common that they were different. The median of the absolute difference between the two bond lengths was 0.35 Å, with a maximum of 0.99 Å. However, in all cases, the hydride ion was directed toward the other Fe ion, in contrast to a terminal binding to the Fe ion in which the hydride ion typically binds trans to the central carbide ion. 33 There was a fair anticorrelation between the absolute difference in the two Fe−H bond lengths and the length of the shorter bond, i.e., that the shortest bonds were found for pairs with  However, the correlation was reduced by structures involving the partial or full formation of H 2 and also by a number of structures in which a third Fe ion also came rather close to the hydride ion (but always with a Fe−H distance larger than 2.02 Å). However, for distances larger than approximately 2.20 Å, such interactions seem to be more occasional than reflecting any specific Fe−H interaction ( Table 2 shows such interactions when they are shorter than 2.3 Å).
The relative energies of the various structures are shown in the last column of Table 2, calculated at the TPSS-D3/def2-SV(P) level of theory. It can be seen that the best structure has one hydride bridging Fe2 and Fe6 on the side toward S3A (Fe2/ 6(3)) and the other hydride bridging Fe3 and Fe7 on the side toward S2A (Fe3/7 (2)). This structure is shown in Figure 4a. As for the structures in Table 1, the hydride ion between Fe3 and Fe7 receives a dihydrogen bond from Arg-96 (the H···H distance is 1.43 Å), whereas the other hydride ion is 3.4 and 3.8 Å from polar protons on Gly-357 and His-195. The two hydride ions are 3.69 Å apart, on different faces on the FeMo cluster, but the Fe3/7 hydride and the S2B are on the same face at a distance of 3.13 Å. The structure is also stabilized by hydrogen bonds between His-195 and S2B (the H···S distance is 2.20 Å) and between S5A and two water molecules (the H···S distances are 2.41 and 2.67 Å).
The structure is quite similar to that suggested by Hoffman and co-workers, 3,28 shown in Figure 2 and included in last two rows in Table 2; all protons and hydride ions are on the same S and Fe ions. However, the proton on S5A points in the opposite direction (viz., toward S3A rather than toward S2B), and the Fe2/6 hydride ion is on the other side of S2B (viz., on the S3A side). In the Hoffman structure, both protons and both hydride ions are on the same face of the FeMo cluster, with H−H distances of 2.04 (the two hydride ions), 2.47, and 2.16 Å. Thereby, the reductive elimination of H 2 can be easily accomplished. 3,28 However, this makes this structure crowded, and it is 27 kJ/mol less stable. This is probably because the hydrogen bond from His-195 to S2B has a worse geometry (the H···S distance is 2.34 Å). The two hydrogen bonds between water and S5A are hardly changed (the H···S distances are 2.41 and 2.69 Å) and neither is the dihydrogen bond between the Fe3/7 hydride and Arg-96 (H···H = 1.44 Å). The Fe2/6 hydride is 3.45 Å from His-195 and 2.79 Å from Arg-96.
Our second-best structure has the hydride ions on the same Fe ions, but the Fe2/6 hydride is on the other side of S2B, directed toward S5A (Fe2/6(5)). On this side, the hydride ion does not interact with any protein residue (it is 2.64 Å from a proton of Arg-96, but this is the same proton that forms a close dihydrogen bond to the other hydride ion, at 1.42 Å distance). This is the best structure from Table 1 (shown in Figure 3) and is thus even closer to the Hoffman structureit differs only in the direction of the proton on S5A. It is 39 kJ/mol less stable than the best structure.
The next structure in terms of energy depends on the BS, but   34 This indicates that H 2 formation is disfavored in these rather crowded structures.
All optimizations of the structures in Table 2 were started from the same BS state. However, it was quite hard to keep the BS state and find the same BS state for all these structures, so several of the structures ended up in differing BS states. The reason for this is primarily because the spin population of some of the Fe ions (typically those involved in the binding of the hydride ions) are low (often below 0.5) and therefore can easily change sign during the optimization of both the wave function and the geometry.
For the best structures in Table 2 and also the Hoffman structure, we performed a thorough study of the best BS state. It turned out to be harder to find the various BS states by swapping the Fe ions 64 with these E 4 structures with two bridging hydride ions than in our former studies. 37,38 Thus, we obtained only 12− 18 of the 35 possible BS states for the four different structures.
Interestingly, this also led to the discovery of a new type of BS state. In all BS states described so far in this study and also in earlier studies, 33,37 four of the Fe ions had a surplus of majority spin (denoted α or positive spin in the following), whereas three had a surplus of minority spin (β or negative). Then, the ground quartet state of the resting E 0 state can formally be obtained by combining one α and three β Fe(III) ions with three α Fe(II) ions and a doublet α Mo(III) ion (−5 − 5 − 5 + 5 + 4 + 4 + 4 + 1 = 3). Likewise, the doublet ground state of E 4 can be obtained by switching the spin of the Mo ion. However, for many binding modes, we found that it was more favorable to have only two Fe ions with β spin.
In fact, the most stable structure in Table 2 turned out to have such a configuration. The best BS state of the Fe2/6(3)−Fe3/ 7(2) structure had a β spin on the Fe1 and Fe4 ions, i.e., the BS-14 state. These two ions also had the largest (absolute) spin populations, both −3.3. Neither of these are involved in the binding of any hydride ion. The third Fe ion that is not involved in hydride binding, Fe5, also had a rather large spin population of 2.5. However, that on Fe3 was slightly larger, 2.6. . This BS state was 19 kJ/mol more stable than the best BS states with three Fe ions with β spins (BS9-126, BS10-147, and BS7-235, which were degenerate; all obtained BS states are shown in Figure S2a). BS-14 was also lowest with the B3LYP functional.
We also tried to find all BS states with only two Fe ions with β spins starting from the BS-14 state (there are 21 7 6 2 = × such states) by systematically swapping the Fe ions. However, most such calculations ended up in states with three β-spin Fe ions, and at the end, we only found two additional states, BS-13 and BS-15 (i.e., both with β spins on Fe1 and on either Fe3 or Fe5). They were 23 and 35 kJ/mol less stable than the BS-14 state.
After the BS investigation, the third best structure had one hydride on Fe3/7(2), whereas the second hydride bridged Fe5 and Fe6. It was also in the BS-14 state, which was 12 kJ/mol more stable than the BS-13 state and 15−17 kJ/mol lower than BS6-156, BS2-234, and BS7-235 (cf. Figure S2c). With B3LYP, Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article instead BS2-234 and BS7-235 were lowest. The Fe5/6 hydride is on the same face as the Fe2/6(3) hydride in the best structure, as can be seen in Figure 4b. The hydride ion is 3.07 Å from the backbone H of Gly-357 but only 2.54 Å from one of the side chain hydrogen atoms of the same residue. The Fe2−H distance is only 2.37 Å, so the structure is not very different from the best structure. The structure was 52 kJ/mol less stable than the best one (both with the BS-14 state). For the Hoffman structure, BS-14 was again found to be lowest in energy. It was 7 kJ/mol more stable than the BS9-137 state, and there were two additional states within 20 kJ/mol, BS7-346 and BS4-356 ( Figure S 2d). However, with B3LYP, instead BS3-124 was lowest, as much as 70 kJ/mol lower than BS-14.
Further Variations of the Best Structures. Both our previous study 33 and this study indicate that the positions of the protons and hydride ions affect each other. Therefore, we performed some additional variations of the best structures. First, we tested to use the other conformation of the proton on S5A, i.e., directed toward S2B rather than S3A, as in the structure suggested by Hoffman and co-workers 3,28 (Figure 2), which was the second best structure in Table 1. The results are collected in the second part of Table 3. Interestingly, this conformation turned out to be more favorable than S5A(3) for the Fe2/6(3)− Fe3/7(2) structure by 8 kJ/mol (still in the BS-14 state). Thus, this new structure, shown in Figure 5a, is the most stable structure so far. This is in contrast to the Hoffman structure, for which it was more favorable to have the second proton on S5A(3), i.e., giving the Fe2/6(5)−Fe3/7(2) structure in Figure  3. The third best Fe3/7(2)−Fe5/6 structure reorganized to the Fe2/6(3)−Fe3/7(2) structure when a proton was moved to the S5A(2) side.
Third, we tested to move the second proton from S5A to S3B, which was also a low-energy structure in Table 1. S3B has also been suggested to be the entrance of protons into the FeMo cluster. 72,73 However, no low-energy state was found for any of the five tested structuresall structures were 86−107 kJ/mol less stable than the best structure (fourth section of Table 3). Moreover, most of the structures reorganized slightly. Fe5/6 became Fe2/6(3) as usual, and in many structures, the Fe3−H bond was cleaved, so that the hydride bound only to Fe7. For the Fe2/6(5)−Fe3/7(2) structure, the proton moved from S3B to Fe6.
Thus, sorting the structures in Table 3 by the TPSS-D3/def2-SV(P) energy (Figure 7), we see that the four best structures all involve hydrides on Fe2/6(3) and Fe3/7(2), representing the four possible combinations of conformations of protons on S2B and S5A, with the S2B(3)−S5A(3) conformations as the best, 5−13 kJ/mol better than the other three. With the larger def2-TZVPD basis set, the order remains, but the variation increases to 6−40 kJ/mol. When evaluated with the B3LYP-D3/def2-SV(P) method (single-point energies), the best structure remains 11−23 kJ/mol more stable than the other structures.
Next on the energy scale comes the first structure with a Fe2/6 hydride on the other side of S2B (Fe2/6(5)). Again, it is most stable with the protons in the S2B(3)−S5A(3) conformations. It is 13 kJ/mol less stable than the best structure (8 kJ/mol with B3LYP). The S2B(3)−S5A(2) conformation is 8−11 kJ/mol less stable with the two functionals. The other two conformations of the two protons are appreciably less stable, 43−53 kJ/mol higher than the S2B(3)−S5A(3) conformations with TPSS (37−45 kJ/mol with B3LYP), and the worst conformation is actually the Hoffman structure.
The lowest structure with the Fe3/7 hydride on the other side of S5A (Fe3/7(3)) is 42 kJ/mol less stable than the best structure, and the structure with both hydrides bridging Fe3 and Fe7 is actually 4 kJ/mol more stable.
We also included the best TPPS structure from our previous investigations 33,34 in the comparison, i.e., with a proton on S2B(3) and three hydride ions on Fe2/6(3), Fe5, and Fe6 (one bridging and two terminal hydride ions). This structure is shown in Figure 6e, and it is called 3H in the following. As in our previous studies, it was found to be most stable in the BS2-234 state (the BS-14 state could not be found, in spite of several attempts). Interestingly, it turned out to be 14 kJ/mol less stable than our best structure with two bridging hydride ions, the S2B(3)−S5A(3)−Fe2/6(3)−Fe3/7(2) structure in Figure 6a, making it the sixth best structure in this investigation. With the B3LYP functional, it was appreciably less stable, 48 kJ/mol above the best structure, reflecting that B3LYP disfavors structures with hydride ions on the Fe ions.
To further check that these energies are reliable, we performed some additional calculations with a larger QM system, including also the side chains of Val-70 and Gln-191, the full side chain of Arg-96, and an extra methyl group for the two His residues (186 atoms in total, shown in Figure S1). The results are also included in Table 3 and Figure 7, and it can be seen that the extended QM system changed the relative energies by up to 22 kJ/mol. In particular, the S2B(3)−S5A(3)−Fe2/ 6(3)−Fe3/7(2) structure remained most stable, but the S2B(3)−S5A(3)-Fe2/6(5)−Fe3/7(2) became second best, only 6 kJ/mol less stable, whereas the 3H structure was destabilized (36 kJ/mol less stable than the best structure) and Hoffman structure remained quite unfavorable (56 kJ/mol less stable than the best structure). All calculations up to now have been performed with the protein outside the 151-atom QM system fixed to the crystal structure (to minimize the local-minima problem, i.e., that different structures end up in different local minima of the Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article surroundings, giving unstable and unreliable energies). On the other hand, there is a risk that this gives a too small flexibility of the calculations (even if the QM system is quite large). Therefore, we have repeated the QM/MM calculations for the best structures allowing the surrounding protein to relax by a MM minimization in each step of the QM/MM geometry optimization. This had a larger effect on the ordering of the structures, although most relative energies did not change by more than 18 kJ/mol. In particular, the S2B(3)−S5A(3)−Fe2/ 6(3)−Fe3/7(2) structure (which was most stable with the fixed surroundings) was destabilized compared to all the other structures, except the S2B(5)−S5A(3)−Fe2/6(3)−Fe3/7 (2) and became the second-best structure. Instead, the S2B(5)− S5A(2)−Fe2/6(3)−Fe3/7(2) structure became the best structure by 13 kJ/mol. Unfortunately, it is hard to see any reason for the stabilization.  ). However, the same residues move also for other E 4 states. In fact, the movement is slightly smaller for the S2B(5)−S5A(2) conformation than for the S2B(3)−S5A(3) conformation, indicating that the former fits the protein structure somewhat better. The latter conformation also gives larger movements of the backbone of Arg-359. Thus, there are some indications that the S2B(5)−S5A(2) conformation may be more favorable in a relaxed surrounding, although it cannot be excluded that this stabilization comes from the local-minima problem, which is much larger when the surroundings are relaxed.  Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article In fact, an even larger stabilization is seen for the Hoffman structure (S2B(5)−S5A(2)−Fe2/6(5)−Fe3/7(2)). With the relaxed surroundings, it is the sixth best structure, 22 kJ/mol worse than the best structure. In this case, the change in the S−H and Fe−H distances are smaller (0.00−0.02 Å), but Val-70 shows an even larger movement (up to 1.4 Å for the HG2 protons), indicating that the surroundings move to relax unfavorable interactions between the protonated S3B and His-195 (the distance between the proton on S2B and the closest HG2 proton on Val-70 is only 1.8 Å; because of the crowding caused by the positioning of all four hydrogen atoms on the same face of the FeMo cluster in the Hoffman structure, the S2B atom and its proton move by 0.76 and 1.07 Å, respectively, compared to the structure with Fe2/6 on the other side of S2B, interfering with the hydrogen bond between His-195 and S2B). Still, the general conclusions of this study remain. With relaxed surroundings, the most stable structure has hydride ions on Fe2/6(3) and Fe3/7 (2), and this structure is 18 kJ/mol more stable than the 3H structure and 22 kJ/mol more stable than the Hoffman structure.
Studies with Other DFT Functionals. So far, we have judged the stability of the structures with the TPSS functional. To investigate how the relative energies depend on the DFT method, we have tested seven additional DFT functionals, viz., the pure generalized gradient approximation (GGA) PBE and M06L and B97D functionals without any Hartree−Fock exchange, as well as the TPSSh, B3LYP, PBE0, and M06 hybrid functionals with 10%, 20%, 25%, and 27% Hartree−Fock exchange, respectively. With these functionals, we compared the energies of the 3H structure with eight of the best E 4 structures found in the present study, all including Fe3/7(2), but differing in the direction of the protons on S2B and S5A, as well as the position of the Fe2/6 bridging hydride.
Thus, we can conclude that for all the pure functionals the most stable E 4 state has two bridging hydride ions. However, for TPSSh, the 3C structure, is approximately 9 kJ/mol more stable, whereas for the other hybrid functionals, the 3C structure is >100 kJ/mol more stable. Thus, after this thorough study, we actually find that with the pure GGA functionals, the energetically most favorable E 4 actually has two bridging hydride ions, in accordance with the ENDOR experiments. 3,29 This solves one important problem in the previous studies based on pure functionals. 34 However, another problem still remains. It is still strongly favorable for H 2 to dissociate from our best E 4 structure, forming an E 2 structure with a proton on S2B and a hydride ion bridging Fe2 and Fe6 (Fe2/6(3)) and an isolated H 2 molecule in a water-like COSMO solvent with a dielectric constant of 80 (the reaction energy is 59 kJ/mol with TPSS-D3/ def2-SV(P)).
Comparison with ENDOR Results. Hoffman and coworkers have studied which DFT structures are compatible with their ENDOR data using an analytical point-dipole model of hydride electron−nuclear dipolar coupling. 30 This analysis is based entirely on the geometry of the two Fe−H−Fe bridges and their relative orientations, employing the fact that the experimental data show purely rhombic dipolar tensors with permuted orientations. They compared five DFT models and showed that the Hoffman model in Figure 2 reproduces the experimental data best, with a mean absolute deviation (MAD) for four angles between the dipolar tensors of only 1°−3°for TPSS and B3LYP structures. The estimated uncertainty in the experimental angles is 5°. 30 We have performed the same analysis 30 on all our structures. For our Hoffman structure, obtained with TPSS, we get an appreciably larger MAD, 9°, mainly because the angle between the two Fe−H−Fe planes (called τ in their article 30 ) is larger, 25°compared to 11°, probably caused by the smaller conformational freedom when the structures are optimized inside the enzyme. However, the ENDOR angles strongly depend on the details of the calculations. The τ angle varies between 3°and 35°when the Hoffman structure is optimized with the eight different DFT functionals. Consequently, the MAD also varies extensively, from 5°with TPSSh to 92°with M06. Likewise, the MAD depends strongly on the BS employed, from 5°to 15°for the 19 obtained BS states, although the variation in τ is quite small, 23°−28°, except for BS9−145, which gives τ = 19°and also the lowest MAD. Moreover, the ENDOR data require that the two Fe ions in one of the Fe−H− Fe bridges have parallel spin, whereas the other two have antiparallel spin. Although this is the case for 16 of the 35 Noodleman BS states, these are not those lowest in energy. The most stable BS state of that type (BS9-137) is the fifth best state, 23 kJ/mol less stable than the best (BS-14 state). The state with the lowest MAD (BS9-145) is 27 kJ/mol less stable than BS-14 and does not have proper spins of the Fe−H−Fe. Hoffman et al. simply selected the lowest BS state with proper alignment, which was within 20 kJ/mol of the lowest state (but they did not consider our best BS-14 state, which is 7 kJ/mol lower than the lowest BS state with three negative spins in our calculations). The MAD with TPSS and BS-14 for both the relaxed protein or the larger QM system is slightly lower, 7°(τ = 23°in both cases).
The MAD is quite insensitive to the position of the two protons. For example, the MAD of the 14 structures in Table 1 differ by only 4°, with the most stable structure (i.e., with the S5A proton pointing toward S3B rather than toward S2B, as in the Hoffman structure) also giving the lowest MAD. The  Table 2 (all with protons on S2B(5) and S5A (3)), the Hoffman-type Fe2/6(5)−Fe3/7(2) structure gives the lowest MAD, 7°with the best BS-14. However, with the same BS used for most of the structures, the MAD is higher, 13°, which is similar to that obtained with several other structures. Fe3/7(3)−Fe2/3 gives MAD = 11°but is high in energy (143 kJ/mol less stable than the best structure). Fe2/6(3)−Fe2/3 gives MAD = 13°with BS-14 but 26°with the standard BS state and is intermediate in energy (87 kJ/mol with BS-14). Likewise, Fe4/5(2)−Fe3/4 and Fe4/5(5)−Fe2/3 give MADs of 13°and 15°but are both high in energy (167 and 187 kJ/mol). However, for the Fe5/6−Fe3/7(2) structure, the MAD is 17°−21°, depending on the BS state, and it is rather low in energy (52 kJ/mol). On the other hand, the energetically best structure in Table 2 (Fe2/6(3)−Fe3/7(2)) gives a poor MAD of 44°(τ = 116°). The same applies to the other Fe2/6(3)− Fe3/7(2) structures in Table 3  Dance 35 has suggested that the Fe spin populations can be compared to ENDOR 57 Fe isotropic hyperfine coupling constants of the E 4 intermediate in the Val70Ile mutant of nitrogenase, although he did not obtain any satisfactory agreement with any of his DFT structures. 74 We have performed a similar comparison of the normalized Mulliken spin populations of all our studied structures (Tables S3−S5) to the normalized hyperfine constants, 35 and we also find no complex that shows agreement with the hyperfine coupling constants. The problem is the same as that encountered by Dance,35 viz., that the five largest spin populations (in absolute terms) are of a more similar size than the measured hyperfine constants. For example, the second largest hyperfine constant is only 77 ± 3% of the largest one, whereas the second largest spin population is on average 91% or 98% of the largest one for the complexes in Tables 2 and 3. In fact, the spin population is below 80% for only six complexes (with 3H as the only low-energy complex). Even worse, the fourth highest hyperfine constant is only 49 ± 3% of the largest one (the third largest one is only indirectly estimated and has a high uncertainty), and only one complex has such a low fourth largest spin population. It is not low in energy and had a large second-largest spin population (96% of the largest one).
Finally, we note that the Mulliken charge on the hydride ions bridging Fe ions (0.01−0.13 e in the various complexes) is only slightly smaller than that on protons bound to S atoms (0.12− 0.24 e), but they always have a significant spin population (0.04−0.09). Thus, in the calculations, there is no large difference between what is called protons and hydride ions (based on formal oxidation-state assignments), as has also been discussed by Dance. 75

■ CONCLUSIONS
We have studied possible atomistic interpretations of the E 4 state in nitrogenase with the aim of finding structures that are in accordance with experimental ENDOR data. 3,29 We have systematically calculated the relative energy of the 66 possible structures with two bridging hydride ions and two protons on sulfides. For the best structures, we have also investigated the preferred sites for the two protons.
Our investigation indicates that the most stable E 4 structures have two hydride ions bridging between Fe2 and Fe6 and between Fe3 and Fe7. These are the same two pairs of Fe ions as suggested by Hoffman and co-workers, based mainly on arguments from experimental observations. 3,28 It is quite satisfying that our calculations confirm that these two pairs actually give the most stable structures also in the QM/MM calculations. This indicates that we are slowly approaching some consensus regarding the mechanism of nitrogenase. However, in the TPSS structure with the lowest energy, the Fe2/6 hydride is on the opposite side of the S2B sulfide, compared to the structure suggested by Hoffman, i.e., directed toward S3A rather than toward S5A (cf. Figures 2 and 5a). In fact, our best structure is 22−56 kJ/mol more stable than the Hoffman structure, a substantial amount.
The energy of the various E 4 structures depend also on the positions of the two protons. Our study indicates that it is most favorable to protonate the S2B and S5A μ 2 sulfide ions. In both cases, the protons can point in two different directions, and these conformations are often quite close in energy (within 20 kJ/ mol), so that the most stable structure depends on details in the calculations (Figure 7). Thus, the S2B(3)−S5A(3)−Fe2/6(3)− Fe3/7(2) structure is normally most stable, but if the surrounding protein is allowed to relax, the S2B (5) The same applies to the Hoffman-type structures with the Fe2/6(5) and Fe3/7(2) hydride ions. All four variants of conformations of protons on S2B and S5A are among the most stable structures. In fact, the original Hoffman structure, involving the S2B(5)−S5A(2) protons, is typically the least stable of the four combinations. Instead, the S2B(3)−S5A(3)− Fe2/6(5)−Fe3/7(2) structure is the most stable of the Hoffman-type structures, except with the large basis set, for which the S2B(3)−S5A(2)−Fe2/6(5)−Fe3/7(2) structure is better. The best of these two structures are always within 6−22 kJ/mol of the best Fe2/6(3)−Fe3/7(2) structure, which is within the uncertainty of the present calculations. Thus, it is possible that the E 4 state may involve the Fe2/6(5)−Fe3/7(2) bridges, which reproduce the experimental angles between the ENDOR dipolar tensors much better than the Fe2/6(3)− Fe3/7(2) structures. However, then, it is likely that the protons are in the S2B(3)−S5A(3) conformations, rather than the S2B(5)−S5A(2) conformations, suggested by Hoffman. 30 Several other studies have addressed the structure of the E 4 state. Recently, Einsle and co-workers have suggested that both hydride ions should bridge the Fe2 and Fe6 ions, after dissociation of the S2B bridging sulfide (without specifying the positions of the two protons). 76 This is rather similar to the first structure (Fe2/6(3)−Fe2/6(5)) in Table 2, which is 99 kJ/ mol less stable than the best structure in that table, although our structure was obtained with S2B still bound to the cluster (and protonated). Hoffman et al. also considered such a structure, but with S2B protonated and binding to Fe6 but not Fe2. 30 It was 16 kJ/mol less stable than their best structure, and it was shown to Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article be incompatible with the ENDOR data (we obtain a MAD of 47°for our Fe2/6(3)−Fe2/6(5) structure). Dance has studied approximately 40 possible structures of the E 4 state. 35 The energetically most stable structures all involved bound H 2 , probably reflecting a tendency of BLYP to overestimate the stability of such complexes. 33,34 The best complex without H 2 (40 kJ/mol less stable than the best H 2 complex) involved three terminal hydride ions, but a complex with one bridging hydride ion and two terminal was only a few kJ/mol less stable. The best complex with two bridging hydride ions (S2B(5)−S3B−Fe2/6(5)−Fe6/7) was 33 kJ/mol less stable, and the S2B(5)−S3B−Fe2/6(5)−Fe3/7(2) structure was 117 kJ/mol less stable. Five additional complexes with two bridging hydride ions (Fe2/6(5) and Fe3/7(2) or Fe6/7) were studied, but all had one or two terminal hydride ions. Thus, there is no overlap with the complexes in this study, although we studied the S2B(3)−S3B−Fe2/6(5)−Fe3/7(2) complex, which differs only in the direction of the proton on S2B, and our results agree that it is high in energy. Dance's results indicate that Fe6/7 is more favorable than Fe3/7, but the reason for this is most likely that he employed a minimal QM-cluster model, which did not include any model of Arg-96, which stabilizes a Fe3/7(2) bridging hydride ion.
Our study shows that for the hydride-bound structures the problem with finding the most stable BS state becomes even more severe than for the less reduced states. The reason for this is that the spin populations on the hydride-bond Fe ions become lower and therefore can more easily change during the optimization of the wave function and of the structure. Moreover, the hydride ions break the approximate C 3 symmetry of the FeMo cluster, so that all 35 states with three β-spin Fe ions need to be considered. Even worse, we find that some states with only two β-spin Fe ions become low in energy, so that also the 21 states of this type should be considered (although in the end, we could only find three of them). It is important to optimize geometries for each BS state, and the different BS states may change relative energies of the various hydride-bound states by up to 50 kJ/mol, even if only low-energy BS states are considered.
Interestingly, our best structure (S2B(3)−S5A(3)−Fe2/ 6(3)−Fe3/7(2)) is 5−48 kJ/mol more stable than our previous best candidate for the E 4 , viz., the 3H structure in Figure 6e with three hydride ions. 33,34 The reason for this discrepancy is mainly the finding of the new BS-14 state, which reduced the energy of our best E 4 structures by approximately 19 kJ/mol. However, it also reflects the failure of our heuristic approach. 33 The present results show that the position of the protons and hydride ions affect each other, whereas in our previous study, we always started from the previous best E n−1 state when finding the best structure of the E n state. Importantly, our new findings solve a problem with DFT studies on nitrogenase. Before, 34 we pointed out that no DFT method gave a most stable structure for E 4 in agreement with the ENDOR experiments suggesting two bridging hydride ions. 3,29 Instead, pure functionals preferred structures with three hydride ions (3H), whereas hybrid functionals preferred structures with a triply protonated carbide ion (3C). 13,26,27,31,32 This made the selection of a DFT method for computational studies of nitrogenase problematic. In fact, Siegbahn and Blomberg has suggested that the resting state needs to be reduced by eight electrons before the actual catalytic E 4 state is reached, for which hybrid functionals suggest a best structure with two bridging hydride ions, in addition to a triply protonated carbide ion, that has moved out of the center of the cluster, and three protonated sulfide ions. 32,77 However, Hoffman et al. has argued that such a structure is inconsistent with experimental data. 30 The present results indicate that pure functionals actually do suggest a structure of E 4 (reduced by four electrons from the resting E 0 state) in accordance with experiments. Pure functionals (and TPSSh) also give geometries of the resting state that are close to the crystal structure. 28, 34 Dance has argued that pure functionals give better energies for reactions related to those involved in nitrogenase. 78 On the other hand, pure functionals still have problems with a too favorable dissociation of H 2 from the nitrogenase models. 34,77 Therefore, further studies are needed before it can be settled which DFT functional gives most reliable results for nitrogenase.
In conclusion, our take-away messages are as follows: • The most stable structures for the E 4 state with pure functionals involve two bridging hydride ions, in accordance with ENDOR data. • The best structures involve bridges between the Fe2/6 and Fe3/7 ion pairs and protons on the S2B and S5A ions.