An assessment of pure, hybrid, meta, and hybrid‐meta GGA density functional theory methods for open‐shell systems: The case of the nonheme iron enzyme 8R–LOX

The performance of a range density functional theory functionals combined in a quantum mechanical (QM)/molecular mechanical (MM) approach was investigated in their ability to reliably provide geometries, electronic distributions, and relative energies of a multicentered open‐shell mechanistic intermediate in the mechanism 8R–Lipoxygenase. With the use of large QM/MM active site chemical models, the smallest average differences in geometries between the catalytically relevant quartet and sextet complexes were obtained with the B3LYP* functional. Moreover, in the case of the relative energies between 4II and 6II, the use of the B3LYP* functional provided a difference of 0.0 kcal mol–1. However, B3LYP± and B3LYP also predicted differences in energies of less than 1 kcal mol–1. In the case of describing the electronic distribution (i.e., spin density), the B3LYP*, B3LYP, or M06‐L functionals appeared to be the most suitable. Overall, the results obtained suggest that for systems with multiple centers having unpaired electrons, the B3LYP* appears most well rounded to provide reliable geometries, electronic structures, and relative energies. © 2012 Wiley Periodicals, Inc.


Introduction
Metalloenzymes, in particular those containing iron, catalyze a broad range of metabolically important reactions within organisms including CAH bond activation, hydrolysis, and DNA repair. [1][2][3][4][5] A detailed understanding of their chemistry and catalytic mechanism provides invaluable insights into their biochemical function as well as those of related species. In addition, it can also provide fundamental chemical insights into, for instance, catalysis, as well as enable the development of new more effective therapeutic treatments. Consequently, they have been the subject of numerous experimental and computational investigations.
For computational investigations on iron-containing metalloenzymes, hybrid density functional theory (DFT) methods, in particular B3LYP, have become the standard tools of choice. [6][7][8][9] It has been widely applied to a variety of such systems and shown to be able to provide useful chemical and mechanistic insights. [4][5][6] This is perhaps surprising given that the B3LYP functional was parameterized based on reference molecules that do not include metals. [10] Indeed, for some metal-containing systems, such approaches have been shown to not give reasonable agreement with experiment or to in fact fail. [11,12] In particular, for Fe-containing systems with near-degenerate states, such methods are often unable to give a qualitatively or quantitatively correct ordering of states. [10] For example, in a study by Ghosh et al. [12] on an FeIII(Por)Cl complex, the B3LYP method predicted that a quartet and not the experimentally observed sextet was the ground state. Unfortunately, accurate determination of the relative energies of near-degenerate states can be essential in elucidating a given systems chemistry. [10,13] The reliability of such methods in describing the relative energies can be sensitive to the amount of exact Hartree-Fock (HF) exchange included. [10] In particular, it has been proposed that for metal-containing systems with unpaired electrons, a HF contribution of 15% provides better accuracy. [10,[14][15][16][17] Indeed, Reiher et al. [10] showed that for B3LYP, the reduction from the widely used value of 20 to 15% (denoted as B3LYP * ) was able to reproduce the energetics of all Fe(II)-S complexes they considered as part of their investigation. In addition, the optimized structural parameters were also in better agreement with experimental X-ray structures. However, more recently, Hughes and Friesner [18] constructed a database of experimental spectra of 57 octahedral first-row transition metal complexes that included V, Ni, Mn, Cr, Fe, and Co. In addition, they also varied the ligands to include examples of MAC, AN, AO, AS, AN, AF and ACl bonds. Then, they examined the effects of reducing the HF contribution in B3LYP from 20 to 15% on its ability to reproduce spin-splitting in such complexes. They concluded that reducing the percent HF contribution did not lead to general satisfactory agreement with experiment. [10] Notably, in all of these systems, any unpaired electrons were formally located on the sole metal-ion centre. That is, the effect of the percent HF contribution on systems containing multiple centers of unpaired electrons has not, to the best of our knowledge, been examined.
Although these differences exist, it is generally accepted that the mechanisms for peroxidation are consistent for all LOXs and are as shown in Scheme 1. The first key step of the process involves a proton-coupled electron transfer from the substrate to the Fe(III)-OH moiety to form a pentadienyl-type radical (I) while the high-spin Fe(III) is reduced to a high-spin Fe(II) centre, as shown by experimental magnetic susceptibility studies, but now ligated to a H 2 O. [28] Subsequently, the pentadienyl radical undergoes an attack by O 2 from the opposite side of the substrate to that of the Fe centre to give the corresponding peroxy-radical (II). Experimental investigations have suggested that this activation of the substrates CAH bond by the Fe(III)-OH group occurs before any kinetically productive interaction between the enzyme-substrate complex and O 2 . [32] Furthermore, the ferric and ferrous centers in I and II, respectively, do not bind O 2 , if at all, before the formation of the pentadienyl radical. [32,33] Indeed, experimental EPR studies have shown that there exists an equilibrium between the pentadienyl and peroxy-radicals, [32,34,35] that is, the O 2 attack is fully reversible. Thus, in the case of LOXs, the enzyme activates the substrate by generating the pentadienyl radical intermediate making it vulnerable to O 2 attack. This contrasts with that observed in most heme-containing enzymes or the family of nonheme aÀketoglutarate-dependent dioxygenases that utilize O 2 to oxidize organic substrates: [2,3] the Fe ion binds the O 2 and facilitates its activation before reaction with the substrate. Thus, in the case of LOXs, the experimental evidence suggests that the key pentadienyl þ O 2 intermediate (II) involves three open-shell centers: the triplet O 2 , the substratederived radical, and the Fe center. Given the high-spin state of the Fe center and the triplet O 2 (i.e., seven unpaired electrons), there are only two possible electronic configurations that allow for spin conservation following the attack of O 2 (i.e., II ! III), illustrated in Scheme 2. Specifically, it involves the arrangement of electrons such that total spins of 3/2 or 5/2 are obtained generating a quartet or sextet state represented as 4 II and 6 II, respectively.
This key intermediate (Scheme 1: II) in the mechanism of the biochemically important LOX family of enzymes presents a clear example of a multicentered open-shell system. The resulting complex likely involves three open-shell centers with a total number of seven unpaired electrons. These open-shell centers being the Fe(II), the pentadienyl radical intermediate (hereafter referred to as AA ), and the molecular oxygen. The catalytically relevant arrangement of these electrons leads to the possibility of having the mechanistic intermediate II with total spins of 3/2 or 5/2; that is, a quartet ( 4 II) or sextet ( 6 II) state (Scheme 2). It is noted that, as for a singlet biradical, a single determinant of KS orbitals (or MOs in the case of HF theory) is incapable of properly describing such multicentered open-shell systems. Instead, an MC-SCF approach must be used. Unfortunately, however, such methods are computationally impractical for the study of most biochemical systems due to the often necessary requirement for large chemical models.
To assess the applicability of commonly used DFT, methods for the study of multicentered open-shell systems, and as a first step toward obtaining an understanding of the chemistry of LOXs, we have examined the performance of a range of hybrid, meta, and hybrid-meta GGA density functionals to reliably determine the structures and energetics of the mechanistic intermediates 4 II and 6 II of 8R-LOX. In addition, we have also considered the effect of varying the %HF contribution in these methods within an ONIOM-type quantum mechanical (QM)/molecular mechanical (MM) model.

Methods
The molecular operating environment (MOE) [36] software package was used to perform all docking and molecular dynamics  (MD) annealing and relaxation of the system. These calculations were done with the AMBER99 force field. [37] Docking The initial crystal structure of 8R-LOX (PDB: 3FG1) was used as a template for docking. All crystallographic waters and counterions were removed. The coordinates of hydrogen were then added using the MOE default method. The substrate, modeled as (4Z,7Z)-1,4,7,10-Undecatetraene, was added in the active site manually and oriented such that the H to be abstracted was within hydrogen distance and interacting with the Fe-OH group ( Figure 2).

MD equilibration
With the substrate docked within the active site, MD simulations were then performed to allow the solvated enzyme-substrate complex to undergo thermal relaxation. In particular, the initial enzyme-substrate complex was solvated with a 7 Å spherical layer of water molecules. To force the system to lie within a volume of space, an ellipsoidal potential wall with a scaling constant of 2 was placed around the solvated enzymesubstrate complex. The damping functional factor included in the electrostatic and van der Waals potentials was set to decay smoothly beyond 8 to 10 Å . Before running the simulations, the geometry of the solvated complex was optimized using the AMBER99 force field until the root mean square gradient of the total energy fell below 0.05 kcal a.u. À1 . The MD simulations were performed under constrained pressure and temperature. The equations of motion were coupled with the Nos e-Poincar e thermostat [38] and the time step for numerical integration was set to 1 fs. Initially, the system was heated from 150 to 300 K for a period of 50 ps, followed by an equilibration period of 100 ps at the constant temperature of 300 K and pressure of 1 atm. A typical structure from the trajectory was then optimized with the AMBER99 force field for the subsequent QM/MM analyses (see below).

QM/MM model
A large active site chemical model was then obtained from the above final optimized structure for all QM/MM calculations. More specifically, it was chosen such that it included the truncated form of the substrate and all active site residues immediately surrounding it, that is, first-shell residues ( Figure 1). In addition, all second-shell residues surrounding the Fe center were retained. It is noted that for added consistency, the same initial structure was used in all optimizations. As a result, any differences observed in the final optimized structures are due solely to the methods and multiplicity chosen.
The a-carbons of each residue was held fixed at their final MM-minimized positions (see above) to ensure integrity of the active site during the calculations. Such an approach has been commonly used in the computational investigation of the catalytic mechanisms of enzymes, and its applicability and reliability have been discussed in detail elsewhere. [39,40] A subset of the complete model centered on the reactive region of the active site was then selected for the high-level QM treatment ( Figure 1: inner circle). Specifically, it consisted of the substrate and those groups directly involved in the reaction. That is, it contained the truncated model substrate, the side chains of His385, His390, His571, and Asn575. In addition, it contained the carboxylate of the terminal Ile694, the attacking O 2 , and the FeAOH center. The O 2 was manually added after the simulation in a cavity consisting of the Gly428 residue proposed to be essential in controlling the O 2 for attack of C8. [41] To form the first intermediate complex (II: Scheme 2), a hydrogen was manually transferred from the substrate to Fe(III)AOH. This complex was the starting point for all further QM/MM calculations.

QM/MM calculations
Combined QM and MM methods in the ONIOM formalism [42][43][44][45][46][47][48][49][50] were performed as implemented in the Gaussian 09 program suite. [51] The optimized structures for 4 II and 6 II (Scheme 2) were obtained using the ONIOM(DFT i /BS1:AMBER) level of theory in the mechanical embedding formulism. [7,8,52,53] DFT i represents the BP86, B3LYP 6 (10% HF contribution), B3LYP * , B3LYP, M06, and M06-L functionals. [10] The combination of basis functions defined by BS1 was the 6-31G(d) basis set on all atoms but Fe, where the LANL2DZþECPs basis set was used. Due to the fixing of atoms within the chemical models, the energies reported are the solely electronic energies. Frequency calculations were performed to validate the nature of the stationary point. Single-point energy calculations on the optimized structures were done at the respective ONIOM(DFT i / BS2//DFT i /BS1:AMBER) level of theory in the electronic embedding formulism. BS2 is defined as the combination of the 6-311G(2df,p) functions plus the LANL2DZ ECPs for the iron. Diffuse functions were not used because, as discussed by Martin et al., [54] their inclusion on Fe is a poor match when used with the underlying 6-31G or 6-311G basis sets. Colour code for residues: included in their entirety (black); modeled as Gly, that is, only the backbone was included with R-groups replaced by a hydrogen (red); modeled as only their R-group, that is, only their C a and side chains included (blue).

Results and discussion
The effect of DFT functional choice on optimized geometries As noted above, the optimized geometries of the quartet and sextet states of the intermediate complex II (i.e. 4 II and 6 II) in the mechanism of LOXs (see Scheme 1) were examined using a variety of density functional methods to describe the high (QM) region of the QM/MM model, the MM method (AMBER) being kept constant. Thus, for simplicity, only the key structural changes observed in the QM chemical model region, which is schematically shown in Figure 2, are discussed herein.
It is noted, however, that the most significant changes on changing the DFT method did occur in the QM region. As can be seen in Figure 2, the metal ion forms a single interaction with each of five amino acid residues; the R-group imidazole of three histidyls, the C-terminus of an isoleucyl, and the Rgroup amide oxygen of an asparaginyl. It also forms a sixth metal-ligand (M-L) interaction with an hydroxyl ion oxygen. Given the distinct nature of the three open-shell centres, a chosen DFT method should give the same or very similar M-L distances for the quartet and sextet states.
As seen in Table 1, however, the MAL bond lengths do differ between the structures of the quartet ( 4 II) and sextet ( 6 II) states. However, the degree of the differences (reported as the absolute value of Dr between complexes) depends on the functional used. For instance, the averages of the absolute differences (Table 1; average |Dr| ) between the optimized bond lengths for 4 II and 6 II as obtained using the BP86, B3LYP 6 , B3LYP * , B3LYP, M06, and M06-L functionals are 0.076, 0.045, 0.009, 0.039, 0.038, and 0.054 Å , respectively. That is, the B3LYP * functional gives a markedly smaller average |Dr| value than the other functionals considered. In contrast, the BP86 and M06-L functionals, that is, those with a 0% HF contribution, have the greatest differences in MAL bond lengths between states. It is noted, however, that M06-L does perform better than BP86 as indicated by it having an average |Dr| value that is $0.2 Å smaller. This is expected given that unlike the latter functional, M06-L was trained with compounds that contain transition metals. [55,56] For B3LYP and M06 (both of which have 20% HF contribution), similar values of average |Dr| were measured of 0.039 and 0.038 Å , respectively. It is interesting to note that although the averages are similar, there are significant differences when comparing individual M-L interactions. With respect to individual interaction distances, all of the B3LYP-based functionals (B3LYP 6 , B3LYP * and B3LYP) predicted the H571 … , N575 … , and I694 … Fe interactions to be shorter, and the H385 … , H390 … , and HO … Fe interactions to be longer in the sextet state compared to that obtained for the quartet state. Similarly, the BP86 functional predicts the N575 … Fe and I694 … Fe distances to be shorter in the sextet state. In contrast, both the M06 and M06-L functionals predict the sextet state to have longer M-L interactions for all ligands, the only exception occuring for the optimized I684-COO -… Fe distance obtained using the M06-L method.
Similar trends were also observed for the differences in the L n ÀMÀL m (where m = n) bond angles (not shown). Specifically, the BP86 method again gave the largest value (2.8 ) for the difference between 4 II and 6 II, whereas the M06 and M06-L functionals gave smaller average values of 1.8 and 2.4 , respectively. The smallest differences were again obtain using the B3LYP-based functionals, B3LYP 6 , B3LYP * , and B3LYP, which gave averages of 1.7, 0.4, and 1.4 , respectively; again B3LYP * gave the smallest error of all functionals considered.
For the triplet O 2 , one of the three ' 'open-shell centres' ' in intermediate II, it was also found that the optimized OAO bond lengths obtained within 4 II and 6 II and the size of the differences between them depended on the choice of functional (Table 1). For instance, for the quartet and sextet complexes, the BP86 functional gave r(OAO) ¼ 1.306 and 1.238 Å , respectively, that is, the O 2 moiety having a longer bond length within the 4 II complex (Table 1). This was also the largest difference (0.068 Å ) observed in lengths between these two states of any functional considered. For the M06 and M06-L functionals, the OAO bond lengths were again predicted to differ significantly though now by 0.047 and 0.057 Å , respectively, with 4 II again having the longer length. In contrast, all of the B3LYP-based functionals had only minor or negligible differences in O 2 bond lengths between the 4 II and 6 II complexes. Indeed, the observed differences for the B3LYP 6 , B3LYP * , and B3LYP functionals were just 0.007, 0.000, and 0.000 Å , respectively. It is also noted that of these three latter  For the substrate radical itself (AA ), modeled as a pentadienyl radical, all DFT methods considered herein gave CAC bond lengths for the quartet and sextet states that were in very close agreement with each other differing by 0.01 Å or less. Similar to that noted above for the L n -M-L m (where m = n) bond angles, the M06 and M06-L functionals gave slightly larger differences in bond angles for AA in 4 II and 6 II for the bond angles obtained using the B3LYP-based functionals (B3LYP 6 , B3LYP * , and B3LYP). More specifically, the use of M06 and M06-L gave differences in the carbon-backbone bond angles of approximately 2.9 and 2.8 , respectively, whereas B3LYP 6 , B3LYP * , and B3LYP all had negligible differences (0.0 ). Much more significant differences among the DFT methods were observed for the carbon-backbone dihedral angles. For instance, the BP86, B3LYP 6 , B3LYP * , and B3LYP methods all gave differences in backbone dihedral angles between 4 II and 6 II that were <2.0 . In contrast, for some of these same backbone dihedral angles, the M06 and M06-L methods gave differences of 14.1 and 12.3 , respectively. This may be due in part to the fact that in the QM/MM optimized structures, it was observed that the active sites of 4 II and 6 II remained fairly consistent, except when using the M06 and M06-L functionals. For the latter two methods, AA shifted considerably in the 4 II complex. As a result, the substrate is able to adopt a markedly more planar geometry than in the corresponding 6 II complexes.
Overall, the above results appear to suggest that within a QM/MM framework, the B3LYP * method is the preferred DFT functional for obtaining consistent optimized structures of a QM region containing multi-open-shell centres.

The effect of DFT method choice on calculated spin densities in 4 II and 6 II
The spin densities (SDs) of various key sites within 4 II and 6 II were then examined. In particular, the SDs on the Fe(II) centre, the oxygen atoms of the O 2 moiety (O1 and O2), and the carbon centre (C8) of AA that is attacked by the O2. In addition, the sum of the SDs on carbons on the pentdienyl radical intermediate itself (AA ) was also considered. All of these values obtained are provided in Table 2.
It can be seen that for both possible multiplicities, the SD on the Fe center is fairly consistent regardless of the choice of DFT method. More specifically, for 4 II and 6 II, it is calculated to lie in the ranges of 3.94-3.80 and 3.87-3.82, respectively. Such SD values for an Fe(II) have been previously concluded to be indicative of a high-spin arrangement. 2 The BP86 and M06-based methods generally gave values near or at the upper ends of these ranges. For example, use the BP86 functional gives values of 3.92 and 3.87 for the Fe(II) centre in 4 II and 6 II, respectively ( Table 2). In contrast, the B3LYP * and B3LYP methods generally gave values toward the lower ends of these ranges.
The calculated SDs on O1 and O2 indicate the presence of the triplet species for all methods used with ranges of 0.92-1.02 and 0.90-1.00, respectively. Notably, the BP86 method gives the lowest SD values for O1 and O2 in both states while the B3LYP 6 and M06-L methods give a low SD value for O2 in the sextet and quartet states, respectively ( Table 2). Meanwhile, the B3LYP * , B3LYP, and M06 methods give SDs on O1 and O2 that are closer to 1 for both 4 II and 6 II.
The most significant variations in calculated SDs, however, were observed in AA itself. As can be seen in Table 2, the B3LYP*, B3LYP, and M06-L methods all give absolute SD values for AA (the sum of the SDs on carbons on the pentdienyl radical) for both 4 II and 6 II that are within 0.02 of 1.0. Furthermore, they also give consistent absolute SD values for the C8 carbon in both states that are approximately 1/3, ranging from 0.34-0.37. In contrast, the BP86 method significantly underestimates the absolute SD on both AA and C8 for the quartet state with values of just 0.55 and 0.20, respectively (Table 2). Meanwhile, the corresponding absolute values in the sextet state using the same functional are 0.97 and 0.35. This difference in calculated SD values is due to the fact that the use of the BP86 method to describe the QM region in the QM/MM model results in greater delocalization of the SD over the amino acid residues ligated to the Fe(II) centre in the quartet system ( 4 II) than in the sextet ( 6 II). This is expected given that DFT tends to suffer from delocalization errors. [57][58][59] In fact, from Table 2, it can be seen that as the %HF is increased, the SD AA in the quartet system approaches a value of 1.00. That is, as the %HF contribution increases, a greater localization of density is observed, which is expected given that HF theory tends to cause overlocalization. [58] Consistency between 4 II and 6 II is reached when the percent HF contribution included is at 15%, that is, for the B3LYP * method. However, it is noted that even with the use of the B3LYP 6 method (i.e., 10% HF contribution), only a modest underestimation of the SD values of AA and C8 in the quartet state compared to the sextet state values is observed (Table 2). More specifically, the absolute SD values of AA and C8 in the quartet state are 0.10 and 0.04 lower, respectively, than obtained for the sextet state using the same functional.
For 4 II, the use of the M06 functional provides an unexpected result. With the use of M06 (which has 20% HF contribution), a marked underestimation of the absolute SDs on AA (0.71) with a concomitant overestimation of the value on C8 (0.47) is observed. This was again found to be due to greater delocalization of SD over the amino acid residues ligated to the Fe(II) center in 4 II compared to that in 6 II. However, as stated above, M06-L (which has 0% HF contribution) gives an absolute SD value for AA for both the quartet and sextet complexes very close to 1.0. Thus, by examining the performance of a range of hybrid, meta, and hybrid-meta GGA density functionals, the above results suggest that at least for ONIOM QM/MM calculations, a reliable description of the SD distribution can be obtained when the B3LYP * , B3LYP, or M06-L functionals are used to describe the QM-region. Thus, it appears that although significant differences in geometry existed between 4 II and 6 II for the B3LYP and M06-L functionals (in comparison to B3LYP*), the proper description of the electronic distribution appears to be less sensitive to geometrical differences, but rather is sensitive to the functional used.
The effect of DFT functional choice on the relative energies of 4

II and 6 II
The reliable and accurate calculation of the thermochemistry of a chemical system is important not only in the elucidation of enzymatic mechanisms but also arguably is a common goal of computational studies. Hence, the ability of the various DFT methods to reliably calculate the relative energy difference between the 4 II and 6 II complexes, that is, between the quartet and sextet states of the 8R-LOX mechanistic intermediate II, was also examined. The results obtained at the ONIOM(DFT i /BS2//DFT i /BS1:AMBER) level of theory within the electronic embedding formulism (EE) are presented in Table 3. That is, relative energies were calculated by performing singlepoint calculations in which the same DFT functional was used to describe the QM-region as for the optimized geometry being used, but now in combination with a significantly larger basis set (BS2; see Computational Methods).
From Table 3, it can be clearly seen that when the B3LYP 6 , B3LYP * , and B3LYP functionals are used to describe the QM region, the energy differences between the quartet and sextet states are within 1 kcal mol -1 with the largest difference (0.6 kcal mol -1 ) occurring for B3LYP. It is noted that both the B3LYP 6 and B3LYP functionals predict the quartet ( 4 II) to lie marginally lower in energy than the sextet ( 6 II). However, as noted in the introduction, Hughes and Friesner 18 found that for several transition metal complexes reduction from 20 to 15%, HF contribution in B3LYP appeared to cause larger errors in the electronic transitions when compared with experimental values. However, on closer inspection of their data, the B3LYP * method in fact appears to give improved results for the Fecontaining compounds. In particular, smaller errors in the e g !t 2g and t 2g !e g spin-forbidden transitions for Fe(TREN-CAM) -3 and Fe(CN) 6 þ3 , respectively, were obtained. In contrast, for the Ni, Cr, and Mn complexes that they considered a reduction to 15%, HF contribution did result in larger errors in the electronic transitions when compared with experimental values.
With the use of the BP86 functional, the quartet is predicted to be 10.9 kcal mol -1 higher in energy than the sextet. Interestingly, with both the M06 and M06-L functionals, energy differences very similar to those obtained with the BP86 functional are observed (see Table 3). In particular, the M06 functional predicts the quartet to be 11.0 kcal mol -1 higher in energy than the sextet, whereas the M06-L functional similarly calculates it to be 10.8 kcal mol -1 higher in energy. Vancoillie et al. [60] have previously shown that for several heme models, the M06 and M06-L functionals overstabilized the high-spin state with respect to the low-spin state. However, as discussed above, the optimized geometries obtained for 4 II and 6 II using the M06 and M06-L functionals to describe the QM-region show significant differences. Thus, it is unclear if the differences in energy between the quartet and sextet states observed for these two functionals are due to the functionals themselves or differences in geometries.
Hence, relative energies were then redetermined via singlepoint calculations on the optimized geometries of the quartet and sextet obtained at the ONIOM(B3LYP*/BS1:AMBER)-ME level of theory; the level that gave the most consistent agreement between the quartet and sextet geometries. That is, relative energies were determined at the ONIOM(DFT i /BS2//B3LYP*/ BS1:AMBER)-EE level of theory and are presented in Table 4.
As can be seen, the energy difference between the two multiplicities has now been reduced significantly to just 0.9 and 1.1 kcal mol -1 for the M06 and M06-L functionals, respectively. It is interesting to note that both functionals predict the quartet state to be slightly higher in energy. This is in contrast to that observed for B3LYP 6 and B3LYP (Table 3). Regardless, it appears that the failure of M06 and M06-L lies predominantly in the determination of the geometry of the quartet system and not preferential stabilization of high-or low-spin. This is further illustrated by the fact that calculation of relative energies at the ONIOM(B3LYP * /BS2//M06/BS1:AMBER)-EE and ONIOM(B3LYP * /BS2//M06-L/BS1:AMBER)-EE levels of theory gives energy differences between the quartet and sextet of 5.2 and 8.2 kcal mol -1 , respectively.

Conclusions
The performance of a range of hybrid, meta and hybrid-meta GGA density functionals to reliably provide geometries, SDs and relative energies of multi-centered open-shell complexes within an ONIOM QM/MM methodology has been examined. More specifically, the ability of the BP86, B3LYP 6 , B3LYP * , B3LYP, M06, and M06-L functionals to provide reliable structures, SDs, and relative energies of a multicentered open-shell mechanistic intermediate complex II in the mechanism of the nonheme iron metalloenzyme 8R-LOX was considered. The latter complex II contains three open-shell centers: an Fe(II), a substrate-derived pentadienyl radical, and a molecular oxygen (O 2 ). From the results obtained, the B3LYP * functional, that is, the B3LYP functional but now with only a 15% rather than a 20% contribution from HF, appears to provide the most reliable geometries. In particular, the use of the B3LYP* method gave the smallest average differences between the catalytically relevant quartet ( 4 II) and sextet ( 6 II) complexes.
In contrast, reliable descriptions of the SD distribution appeared to be obtained using the B3LYP * , B3LYP, or M06-L functionals. Thus, while the B3LYP and M06-L functionals appeared less suited for the proper description of the geometries of 4 II and 6 II, they are capable of properly describing their electronic distribution.
In the case of the relative energies between 4 II and 6 II, the use of the B3LYP * functional at the ONIOM(DFT i /BS2//DFT i / BS1:AMBER) level of theory within the electronic embedding formulism provided a difference of 0.0 kcal mol -1 between the two states. However, B3LYP 6 and B3LYP also predicted differences in energies of less than 1 kcal mol -1 . The erroneously large differences in relative energies predicted using the M06 and M06-L functionals, when using geometries obtained using the same functionals, were found to be due to errors in their underlying optimized geometries. The use of more reliable structures of 4 II and 6 II lead to M06 and M06-L predicting only small relative energy differences between the two complexes.
Overall, the results obtained suggest that for systems with multiple centers having unpaired electrons the B3LYP * appears most well rounded to provide reliable geometries, electronic structures, and relative energies.