1. Introduction

Modeling the reactivities of biological systems is currently one of the long-range goals of computational chemistry. In the field of bioinorganic chemistry, which includes metalloenzyme chemistry, because of improved spectroscopic techniques, the elucidation of the structure of active sites of many enzymes and the nature of the key intermediates in enzyme-catalyzed processes, and an ever-increasing crystallographic database of protein structures, sufficient material has been accumulated to allow the study of many enzymatic catalytic processes in detail at the molecular level [1].

In the last few decades, quantum-mechanical calculations of many chemical reactions including organic, inorganic, and organometallic reactions have provided new qualitative and quantitative insights that have not been available from experimental studies alone. For instance, it was rare ten years ago that organometallic chemists envisioned the structure and energetics of transition states in proposing a reaction mechanism. Today, they often visualize the transition-state structure and collaborate with quantum chemists to find such transition states computationally in order to substantiate a proposed mechanism. The Nobel Prize in Chemistry awarded to John A. Pople and Walter Kohn in 1998 represents the recognition that electronic structure calculation is making a unique and essential contribution to chemistry.

Computational modeling of the structure and functions of metalloenzymes has been slow, because 1) not enough details of relevant reaction mechanisms have been known experimentally for metalloenzyme systems, 2) the electronic structure of metal-containing systems has not been clearly understood, and 3) computational requirements for the modeling of metal-containing systems are very large. However, with problems 1 and 2 being resolved as mentioned above, and computers and codes becoming faster, this field has gained momentum in the last few years. Using rather simple models of enzymatic systems, quantum-mechanical calculations have already led to new insights regarding the behavior of electrons during metalloenzyme reaction processes. In the next section of this review, we restrict our consideration to one case study we performed, the reactivity of methane monooxygenase (MMO) and ribonucleotide reductase (RNR), and discuss some of our findings from the model calculations.

However, the reality of such modeling is still very much restricted by practical computational limitations. As indicated in the next section, computational approaches require the use of relatively simple models of the complicated enzyme system. This is because the electronic structure calculation required to describe such reaction systems demands a great deal of computer time, and our present computers are not fast enough to allow the use of larger and more realistic models. However, as is discussed in Section 3, it appears that even if the computers were faster, one could not apply the theoretical method used in the current model calculations to larger and more realistic models. To achieve realistic results, a highly accurate method is certainly required for the part of the molecule where the reaction is taking place, but it may be possible to adopt less accurate and therefore less computer-demanding methods for the part of the molecule not directly involved in the reaction. We have recently implemented a hybrid approach in which different methods are combined for computations pertaining to one large molecule or a molecular system. In Section 3, we review this “onion-like” method, which we named ONIOM, which facilitates such computations.

Although the ONIOM method has not yet been used for the study of metalloenzyme systems, in Section 4 we briefly discuss its applicability to more realistic models of enzymatic reactions.

2. Density-functional studies of the reactivities of methane monooxygenase and ribonucleotide reductase

Binuclear Fe-containing metalloenzymes, such as hemerythin [2], ribonucleotide reductase (RNR) [3], methane monooxygenase (MMO) [3(a, b), 4], ⁹-desaturase [5], phenol hydroxylase [6], xylene monooxygenase [7], alkane hydroxylase [7], toluene 2-monooxygenase [8], toluene 4-monooxygenase [9], purple acid phosphatase [10], rubrerythrin [11], and ferritin [12], have been the focus of intensive research for the past decade or so. Here we discuss the geometric and electronic properties of the active sites of MMO and RNR, the plausible mechanism of dioxygen coordination and activation on the MMO and RNR, and the mechanism of hydrocarbon hydroxylation by MMO.

Methane monooxygenase (MMO) is an enzyme that catalyzes the methane oxidation reaction, i.e., conversion of the inert methane molecule to methanol [1]:

CH4 + O2 + NAD†(P)H + H+ CH3OH + NAD(P)+ + H2O.

(1)

^† Nicotinamide adenine dinucleotide

The best-characterized forms of the soluble MMO (sMMO) [1, 13] contain three protein components: hydroxylase (MMOH), the so-called B component (MMOB), and reductase (MMOR), each of which is required for efficient substrate hydroxylation coupled to NADH oxidation. The hydroxylase, MMOH, which binds O₂ and the substrate and catalyzes oxidation, is a dimer, each half of which contains three subunit types (, ß, ) and a hydroxy-bridged binuclear iron cluster. In the resting state of MMOH (MMOH_ox), the diiron cluster is in the diferric state [Fe^III–Fe^III], and can accept one or two electrons to generate the mixed-valence [Fe^III–Fe^II] or diferrous [Fe^II–Fe^II] states, respectively. The diferrous state of hydroxylase (MMOH_red) is the only one capable of reacting with dioxygen and initiating the catalytic cycle.

RNR catalyzes the reduction of ribonucleotides and generates deoxyribonucleotides [14]. In the literature there exist at least four classes of RNR, each of which has its unique composition and cofactor requirement [15]. The best-characterized class (I) has an ₂ß₂ structure with ₂ and ß₂ dimers, referred to as R1 and R2 units. Each R2 unit contains a binuclear non-heme Fe active site. This class of RNR requires the presence of O₂ to generate the tyrosyl radical (Tyr) and binuclear active site via oxidation of a binuclear ferrous center of R2 to an oxo-bridged binuclear ferric site:

R2–apo‡ + Tyr122–OH + 2FeII + e– + O2 + H+ R2–(FeIII–O–FeIII) + Tyr122–O• + H2O.

(2)

^‡ Apo-enzyme

Extensive studies of the core structures of the active sites of MMOH and R2 show some sequence homology between these two enzymes, and a similar flexibility of the ligand environment of the Fe-centers [1, 16]. This ligand flexibility has been postulated to be one of the important factors for the proper functioning of the enzymes [17]. Indeed, X-ray studies of the core structure of oxidized MMOH (MMOH_ox) isolated from M. capsulatus show that at 4°C the two Fe atoms are triply bridged by one hydroxo and two µ-1,2-carboxylate ligands with an Fe–Fe distance of ~3.4 Å [18]. Each Fe-center has one histidine ligand. In addition, Fe¹ has one terminal aquo and one carboxylate ligand, while Fe² has two carboxylate ligands, as shown in Scheme 1. On the other hand, the structure recorded at –160°C shows an aquo bridge replacing one of the carboxylate bridges; the resulting structure has one µ-OH, one µ-H₂O, and one carboxylate bridge with a shorter Fe–Fe distance (~3.1 Å) [18(b)]. Similarly, the oxidized binuclear active site of the R2 subunit of RNR from E. coli (R2_met) shows one µ-oxo and one 1,2-carboxylate bridge with an Fe–Fe distance of ~3.2 Å, and one histidine ligand for each Fe-center [19]. Fe¹ has a terminal aquo and a chelating carboxylate ligand, whereas Fe² has two monodentate terminal carboxylate ligands (as in MMOH), and one terminal aquo ligand.

Scheme 1

As shown in the literature [1], the oxidized forms of MMOH and R2 (MMOH_ox and R2_met), including two ferric atoms, Fe^III, are the resting state of these enzymes. Only their two-electron reduced forms, MMOH_red and R2_red, with two ferrous Fe^II iron centers are capable of reacting with O₂—the reaction that initiates both of the catalytic cycles described above. Structural studies [1, 18–20] demonstrated that two-electron reduction of MMOH_ox and R2_met dramatically changes the ligand environment of the Fe centers. Indeed, for MMOH, during reduction two hydroxo/aquo bridging ligands move out, and one of the carboxylate ligands of Fe², Glu243, shifts to form a monodentate bridge between the two metals and coordinates with the Fe²-center in a bidentate manner [18(c)]. Therefore, the two Fe atoms are changed from six-coordinate (6C) to five-coordinate (5C), with one vacant site for each Fe. Spectroscopic studies are in accord with such a (5C, 5C) assignment for the reduced state of MMOH [21]. Similarly, in R2_met the two-electron reduction leads to dissociation of two oxo/aquo bridging ligands [20]. However, upon reduction of R2_met, the terminal ligand Asp84 shifts from chelating to monodentate terminal bridging to Fe¹, a shift which has no analog in the counterpart reduction of MMOH. The carboxylate ligand Glu238 in R2, which moves from terminal monodentate in Fe² to a bridging position, also gives a bidentate bridge rather than monodentate, as in MMOH. Thus, according to crystallographic studies [20], R2_red has two approximately equivalent 4C Fe-centers. However, it is believed that one of the Fe-centers, Fe¹, also has an aquo ligand (see Scheme 1). This is consistent with 1) spectroscopic data [22] indicating that Fe¹ and Fe² atoms in R2_red are 5C and 4C, respectively, and 2) the reduced binuclear active site of manganese-substituted E. coli, which has been found to have a water ligand bound to Mn [1, 23].

Thus, the data described above clearly show that the ubiquitous carboxylate ligands (glutamates and aspartates) of MMOH and R2 may play a special role in the proper functioning of these enzymes. Indeed, the flexibility of these ligands allows them, depending on enzymatic needs, to coordinate with the Fe-centers as bidentate ligands (bridging or chelating) when saturation of the first coordination sphere of Fe-centers is required, and as monodentate (terminal) ligands when one or more vacant coordination sites in the Fe-centers are needed to enable a certain reaction step to take place. Therefore, modeling of the flexibility of the MMOH and R2 ligand environment and the factors affecting this flexibility are important tasks. Modeling the flexibility is the first task of our computational studies.

It has been established in the literature [1] that MMOH_red reacts very rapidly with O₂, forming a metastable, so-called compound O that converts spontaneously to another compound designated as P, (Scheme 2). Spectroscopic studies [24] indicate that compound P is a peroxide species, in which both oxygens are bound symmetrically to the iron atoms. The structure of this compound is unknown. Compound P converts spontaneously to compound Q, which has been proposed to contain two antiferromagnetically coupled high-spin Fe^IV centers. EXAFS and spectroscopic studies [25] of compound Q have demonstrated that it has a diamond-core, (Fe^IV)₂(µ-O)₂ structure with one short (1.77 Å) and one long (2.05 Å) Fe–O bond per Fe atom and a short Fe–Fe distance of 2.46 Å. Compound Q has been proposed to be the key oxidizing species for MMO [1].

Scheme2

In contrast with MMOH, in R2 the second intermediate, X, of the two intermediates U and X observed for the reaction of O₂ with R2, is found to be a spin-coupled [Fe^IIIFe^IV] complex containing two or three µ-oxo-bridges, but no peroxo-level intermediate has been observed. The intermediate U is believed to be a protonated tryptophan radical [1(b)]. Therefore, the second task of our computational studies is to model the mechanism of the dioxygen coordination and O–O bond activation on MMOH and R2 enzymes.

The third important problem related to MMO and RNR is to elucidate the mechanism of the oxygen-atom incorporation into the substrate from the key species of the reaction of MMOH_red and R2_red with O₂, to form compounds Q and X, respectively. In the literature, several mechanisms [1, 26] have been proposed for the reaction of intermediate Q with alkanes; these mechanisms can be divided into two different classes: radical and nonradical. The radical mechanism begins with abstraction of the hydrogen atom from the substrate to form QH (hydroxyl-bridged Q-compound) and a free alkyl radical, whereas the nonradical mechanism implies a concerted pathway, occurring via a four-center transition state and leading to the “hydrido-alkyl-Q” compound. The latest experimental studies [27] show that the mechanism of this reaction is more complex, and could even be different for the enzyme isolated from different organisms. Radical clock substrate probe studies with sMMO isolated from M. capsulatus (Bath) afforded no evidence for the formation of a substrate radical, whereas for the enzyme isolated from M. trichosporium OB3b, evidence for a substrate radical was detected. Thus, the elucidation of the validity of different mechanisms including the structure of assumed intermediates and transition states is an open question and requires extensive further studies.

To provide some insight into these extremely important problems, we have carried out comprehensive computational modeling studies. However, in order to perform high-level quantum chemical calculations on these systems, several nontrivial questions must first be answered.

The first question concerns the choice of a reasonable model for the enzymes incorporating all available experimental findings. For example, according to the experimental results and previous modeling work [28–30], a reasonable model of MMOH_red and R2_red should include two Fe–Fe bridging carboxylate ligands, and one imidazole and one monodentate-coordinated carboxylate ligand for each Fe center. Furthermore, imidazole rings of the His ligands should be located cis to each other. Here, we use two models, shown in Scheme 3. The “small” model includes two Fe–Fe bridging C₁ carboxylates, one NH₂, and H₂O ligands coordinated with each Fe, which respectively model His and monodentate carboxylate residues. The “medium” model used here replaces NH₂ ligands with an imidazole ring, and H₂O ligands with the C₁ carboxylate ligands in the small model. This model additionally includes a water molecule coordinated with one of the Fe (Fe¹) centers.

Scheme 3

Today, these smaller model compounds, in which O–O bond cleavage and oxygen atom incorporation into the substrate take place, can be studied using a high-level computational method. However, one should note that neither model takes into account the protein environment, and, consequently, the existence of different amino acids on the active sites of MMOH and RNR. In other words, in these models, MMOH_red and R2_red species are different from each other ONLY in the coordination mode of one of the bridging carboxylate ligands; in MMOH_red both bridging carboxylate ligands are µ-1,2-type bridges between two Fe centers, while in R2_red one of these carboxylates is a µ-1,1-type bridge. Henceforth, therefore, we call the structures with two µ-1,2-type carboxylate bridges “R2-like” structures, while complexes with one µ-1,2-type and one µ-1,1-type bridge are labeled as “MMOH-like” complexes.

The second important question in these studies concerns the spin state of the system. For example, spectroscopic studies [4] indicate that the binuclear iron cluster of MMOH_ox (diferric) has two S = 5/2 iron atoms antiferromagnetically spin-coupled (AFSC) to yield a diamagnetic center with the exchange coupling constant J between two Fe-centers of 7 ± 3 cm^–1 [4]. However, diferrous MMOH (MMOH_red) exhibits ferromagnetic spin coupling (FSC) of two S = 2 ferrous iron atoms to give an S = 4 ground state [1, 21]. Thus, MMOH_ox and many other intermediates (see below) are AFSC diamagnetic species, while MMOH_red is an FSC paramagnetic species. In other words, during the hydrocarbon hydroxylation reaction by MMO, one may expect spin-crossing of the system and switching between antiferromagnetic and ferromagnetic couplings among the expected intermediates and transition states. The spectroscopic picture for R2_red is also complex. On the basis of Mössbauer studies, the two Fe-centers in R2_red are high-spin ferrous ions [31]. EPR studies of R2_red show a very weak integer spin signal, considered to derive from a small fraction of molecules having FSC sites [32]. MCD studies show a paramagnetic center with a saturation behavior indicating two Fe-centers with M_S = ±2 at the ground state. However, a spin Hamiltonian analysis of the saturation magnetization behavior indicates that the two Fe atoms are weakly (J ~ 0.5 cm^–1) AFSC [33].

The third question involves the choice of adequate computational methods and approximations. Here, the problem should be divided into two parts. The first is an intrinsic problem of how to represent the open-shell low- spin coupling (2S + 1 = 0) of the two paramagnetic centers using the single-determinant methods (such as DFT); the “restricted” calculations for the low-spin state give the closed-shell singlet state, which cannot represent an open-shell atom. An “unrestricted” calculation with low M_S value gives a heavily spin-contaminated state, or the total wave function becomes symmetry-broken and does not represent a true eigenstate. Here one should use a multireference method (such as MCSCF or CASPT2), which, unfortunately, is impractical for such large systems. Therefore, it is sometimes more practical (when the magnitude of the spin coupling between the two centers is not strong) to ignore the antiferromagnetic nature of these systems, and to perform spin-unrestricted open-shell single-determinant calculations for ferromagnetically coupled high-spin states. This approach also retains the proper spins on individual metal (Fe) atoms. Since the magnitude of the spin coupling between the two centers of MMOH and R2 is very small, as discussed above, we expect that the mechanisms of the reactions studied below are not overly influenced by the antiferromagnetic nature of the complexes.

Using the answers to the three questions as described above, we have carried out the model calculations and describe the results as follows. In the first subsection which follows, we discuss the mechanisms and factors affecting the flexibility of the ligand environment of MMOH_red and R2_red. In the next subsection, we discuss the mechanism of dioxygen coordination and activation on MMOH_red and R2_red. The third subsection is devoted to studies of the reaction mechanism of compound Q with the CH₄ molecule, leading to the methanol complex. The last subsection presents a perspective for future studies.

Flexibility of the ligand environment of MMOH and R2
Several attempts have been made to elucidate the factors affecting ligand flexibility and the reactivity of MMOH and R2. Recently, Nordlund and coworkers [17], using mutation experiments, demonstrated that substrate (azide) coordination with R2_red could cause a 1,2-carboxylate shift of one of the bridging carboxylate ligands, Glu238. These authors believe that such a carboxylate shift (of Glu238) could be a key feature in understanding the dioxygen activation mechanism at the Fe-centers. Also, extensive biomimetic studies of carboxylate coordination modes in diiron complexes, as well as O₂ coordination with these biomimetic complexes by Lippard and coworkers [34] have demonstrated that 1) having bulky steric ligands [such as XDK, m-xylylenediamine bis(propyl Kemp's triacid)imide, derivatives] in the bridging carboxylate ligands could facilitate a µ-1,1 µ-1,2-carboxylate shift; and 2) the 1,2-carboxylate shift could be a rate-determining step for the entire substrate (O₂) coordination and ligand rearrangement process [34].

In connection with these experiments, we have also attempted to investigate the carboxylate shift between two Fe centers (called the “1,2-carboxylate shift”) and the bidentate monodentate carboxylate rearrangement within one metal center [29]. In these studies we have used the B3LYP density functional theory (DFT) method [35], which is discussed more fully on urls 381 and 382, and have chosen the “medium” model of MMOH_red and R2_red, presented in Scheme 3. Since the exchange coupling constant of R2_red is extremely small, we have considered both MMOH_red and R2_red as ferromagnetically coupled high-spin species, and we have studied them at their 2M_S + 1 = 9 spin states.

Geometries of the calculated equilibrium structures and transition states (TSs) are shown in Figures 1 and 2, respectively. In Figure 3, we schematically present the calculated mechanisms of the µ-1,1 µ-1,2 carboxylate shift, and the bidentate monodentate carboxylate rearrangement in MMOH-like and R2-like complexes, together with their corresponding relative energies.

Figure 1 Figure 2 Figure 3

As seen from Figure 1, we have found two sets of minima [29]: a) structures 1 and 3, both with one monodentate (µ-1,1) and one bidentate (µ-1,2) Fe–Fe bridging carboxylate ligand (MMOH-like structures); and b) structures 2 and 4, both with two bidentate (µ-1,2) Fe–Fe bridging carboxylate ligands (R2-like structures). Structures 1 and 3, as well as 2 and 4, differ from each other only in the coordination mode of the terminal carboxylate ligand (which models Glu114 in MMOH_red and Asp84 in R2_red) located on Fe¹; the terminal carboxylate ligand binds to the Fe-center monodentately in structures 1 and 2, but bidentately in structures 3 and 4. Thus, horizontally (in Figure 3), processes 1 2 and 3 4 involve migration of one of the carboxylate ligands between the two Fe centers, from being a monodentate bridge (µ-1,1) between the two metals (as well as coordinating with the Fe²-center in a bidentate manner, left), to form a bidentate (µ-1,2) bridge between the two irons (right). This is the so-called 1,2-carboxylate shift. During the first process (1 2), the terminal carboxylate of Fe¹ is monodentately coordinated with the metal center, while during the second process (3 4) this carboxylate coordinates bidentately (chelating) with the Fe¹-center. In other words, during the first process (1 2) the coordination numbers of the Fe centers change from (5C, 5C) in 1 to (5C, 4C) in 2, while during the second process (3 4) the coordination numbers of the Fe centers change from (6C, 5C) in 3 to (6C, 4C) in 4. Meantime, vertically in Figure 3, the processes 1 3 and 2 4 correspond to a monodentate (top) bidentate (bottom) rearrangement of the terminal carboxylate ligand of the Fe¹-center. During these processes, the coordination number of Fe¹ increases by one, from 5C in 1 (and 2) to 6C in 3 (and 4).

First, we discuss processes 1 2 and 3 4, corresponding to the 1,2-carboxylate shift. Comparison of the energies of the 1 and 2 pairs and of the 3 and 4 pairs (Figure 3) shows that R2-like structures 2 and 4, with lower coordination numbers, are energetically preferred over their MMOH-like analogs, 1 and 3; structure 2 is calculated to be 7.9 kcal/mol more stable than structure 1. Similarly, structure 4 is 4.9 kcal/mol more stable than structure 3. The obtained difference in the energies of processes 1 2 and 3 4 can also be due to the ligand environment of the Fe centers; indeed, in the first process we go from (5C, 5C) in 1 to (5C, 4C) in 2, whereas in the second process we go from (6C, 5C) in 3 to (6C, 4C) in 4. Starting from the most stable (2) with the lowest (5C, 4C) coordination, one extra coordination on Fe¹ gives 1.2 kcal/mol of destabilization for 4 with (6C, 4C), one extra coordination on Fe² gives 7.9 kcal/mol of destabilization for 1 with (5C, 5C), and two simultaneous extra coordinations give 6.1 kcal/mol for 3 with (6C, 5C), which is lower than the sum of the two single coordinations (1.2 + 7.9) by 3.0 kcal/mol due to a cooperative effect.

These changes in the ligand environment also have an effect on the metal–metal distance. From MMOH-like structures 1 and 3 (left) to R2-like structures 2 and 4 (right), the Fe¹–Fe² distance increases by ~0.5 Å because of the loss of one (carboxylate) bridge. This is in excellent agreement with the difference in the Fe–Fe distance for the crystallographically determined structures of MMOH_red and R2_red [18(c), 20] and also with recent mutagenic studies involving µ-1,1 µ-1,2 carboxylate shifts [16, 17].

Structures 1 and 2, as well as 3 and 4, are separated from each other by small energy barriers, 0.7 and 1.6 kcal/mol, respectively. The transition states corresponding to these barrier heights, TS12 and TS34 (Figure 2), each have only one imaginary frequency, 88i and 84i cm^–1, respectively, which, according to our normal mode analysis, correspond primarily to the 1,2-carboxylate shift. The first process (1 2), during which Fe¹ remains 5C while Fe² changes its ligand environment from 5C to 4C, is kinetically and thermodynamically more favorable than the second process (3 4), during which Fe¹ remains 6C but Fe² changes its ligand environment again from 5C to 4C. Our calculations clearly show that rearrangements 1 2 and 3 4 are exothermic processes and could take place with very small energy barriers. The key point of our calculations is that the 1,2-carboxylate shift involved in MMOH_red and R2_red structures is an easy process both thermodynamically and kinetically. This is consistent with the experimental observations described above [18–22].

Next, let us discuss the monodentate bidentate rearrangement of the terminal carboxylate ligand in Fe¹, i.e., processes 1 3 and 2 4 for MMOH-like (left) and R2-like (right), respectively. As one can see from Figure 3, for an MMOH-like structure, the carboxylate ligand tends to be chelating (structure 3) rather than monodentate terminal (structure 1); the process 1 3 is calculated to be slightly exothermic by 1.8 kcal/mol, and occurs with a small (1.0 kcal/mol) energy barrier at transition state TS13, calculated from 1. We expect that small variations in the reaction conditions (pH, temperature) may easily reverse the preference observed so far. The computed energy difference (1 kcal/mol) is rather small, and the important fact is that the interconversion is nearly thermodynamically neutral and easily takes place kinetically because of the small barrier.

On the contrary, in the case of R2-like structures, the carboxylate ligand modeling Asp84 tends to be monodentate terminal (2) rather than chelating (4) to the Fe¹-center. The process 2 4 is found to be slightly endothermic by 1.2 kcal/mol, and takes place with a 1.4-kcal/mol barrier at the transition state TS24 (Figure 3). The present result agrees well with experimental observations [6] showing that in the reduced form of R2, the terminal ligand Asp84 tends to be monodentately coordinated with the Fe¹-center.

From these results and discussions, we may conclude that the MMOH-like (1 and 3) and R2-like (2 and 4) structures are energetically very close to each other and are separated by only small energy barriers [29]. Structures 2 and 4 (corresponding to R2_red) with lower coordination numbers in the Fe²-center are found to be a few kcal/mol more stable than those with higher coordination numbers (1 and 3, respectively). Thus, both the 1,2-carboxylate shift (i.e., the shift of one of the bridging carboxylate ligands from µ-1,1 to µ-1,2 bridging mode between the two Fe centers) and the monodentate bidentate rearrangement of terminal carboxylate ligands take place very easily in these systems. Moreover, these reactions can take place reversibly under proper experimental conditions. Our results are consistent with the available experimental data [18–22] showing high flexibility of the carboxylate ligands, as well as a nonrigid core structure, in MMOH_red and R2_red.

Mechanism of dioxygen coordination and O–O bond activation on MMOH_red and R2_red
Dioxygen coordination and O–O bond activation on MMOH_red and R2_red are important steps in the conversion of methane to methanol and the formation of a stable Tyr radical, respectively. This is a nontrivial, multistep process. X-ray studies have revealed [1] that the dioxygen molecule most likely attacks the enzyme from the O-ligand side (opposite the protein residues His147 and His246) because of the existence of the substrate coordination pocket. In the literature, four different coordination positions for the O₂ molecule were discussed (see Scheme 4) [18, 36, 37]. As discussed above, MMOH_red reacts very rapidly with O₂ and forms compound O, which is spontaneously converted to compound P. Compound P is spontaneously converted to compound Q with two antiferromagnetically coupled Fe^IV centers. While spectroscopic studies [24] indicate that compound P is a peroxide species in which both oxygen atoms are bound symmetrically to the iron atoms, the specific coordination mode of O₂ still remains disputable.

Scheme 4

On the basis of available experimental data, as well as our results presented above, which indicate that a) MMOH-like and R2-like complexes differ ONLY in the coordination mode of one of two bridging carboxylate ligands, and b) R2-like structures are a few kcal/mol more favorable than MMOH-like structures, and those can be converted to one another with a small energy barrier, we have proposed (see Schemes 5 and 6) the mechanisms for dioxygen coordination to MMOH-like and R2-like complexes, as well as homolytic and heterolytic O–O bond activation in these compounds.

Scheme 5 Scheme 6

As seen in Scheme 5, the dioxygen molecule can coordinate with either MMOH-like or R2-like complexes I to give O-type complexes (II_MMOH and II_R2, respectively). One may expect that these complexes, like their respective starting complexes I_MMOH and I_R2, are energetically close, separated by a small barrier, and can easily rearrange to one another and into peroxo-type complexes (III_MMOH and III_R2, respectively), in which the O₂ molecule is symmetrically coordinated with two Fe-centers in the cis-µ-1,2 manner. Structures III are, most likely, the intermediate P proposed by experimentalists. From these complexes, the reaction may proceed via the two different mechanisms (homolytic and heterolytic) presented in Schemes 5 and 6, respectively.

The first step in the homolytic O–O activation mechanism is O–O bond cleavage, followed by the insertion of O-atoms between the two Fe centers to form Fe(µ-O)₂Fe-type bis-µ-oxo-complexes V. Our studies indicate that this process occurs with a large activation barrier [38] if during the reaction all Fe–O(bridging carboxylate) bonds remain intact. However, the existing terminal water ligand on the Fe¹ center facilitates the opening of one of the “legs” of the bridging carboxylate ligand on III_R2, and leads to the formation of complex IV. Similarly, structure IV can easily be generated from intermediate III_MMOH by transferring the µ-1,1-carboxylate from the Fe¹(µ-O)Fe² position to the terminal water molecule. Later, intermediate IV rearranges to complex V with a small O–O bond activation barrier. In its turn, complex V rearranges to the Q-type structures VI_MMOH and VI_R2 by transferring one of the carboxylate legs from the terminal water ligand to a µ²-O position between two Fe-centers with a weak interaction with the Fe¹ center.

Support for the heterolytic O–O bond cleavage (see Scheme 6) comes primarily from experimental studies [39] which show that two protons are required in the O₂ reaction, one for the formation of the P-type intermediates (III_MMOH and III_R2) and the other for the conversion of the P-type intermediate to the Q-type intermediates. Lipscomb and co-workers [1, 40] have proposed that a water molecule dissociates as the O–O bond of the P-type intermediate is heterolytically cleaved in the formation of the Q-type structures.

In order to study in detail the mechanisms proposed above, one must formulate the expected spin states of the compounds O, P, and other (not observed) intermediates. The ground state of the O₂ molecule is a triplet ³ state with two unpaired electrons. In addition, the ground electronic state of the other reactants, complexes I_MMOH and I_R2, has eight unpaired electrons each from the two ferromagnetically coupled diferrous Fe-centers. Therefore, one may expect that the O-type compounds II_MMOH and II_R2 are a result of the interaction of triplet O₂ with nonet I_MMOH and I_R2. The resultant complex, [Fe^II(O₂)Fe^II], may have various electronic states, among which the ⁷A, ⁹A, and ¹¹A states should be energetically most favorable.

As expected, the interaction of the triplet O₂ molecule with nonet I_MMOH and I_R2 is extremely weak, and most likely leads to a very shallow minimum corresponding to an ¹¹A state of [Fe^II(O₂)Fe^II]. The latter could convert almost spontaneously to either mixed-valence species [Fe^II(O₂^–)Fe^III] or the compound [Fe^III(O₂^2–)Fe^III], respectively, via extremely fast one- or two-electron transfer from the Fe-centers (d-orbitals) of II_MMOH and II_R2 to the O₂ molecule (* orbitals). Several electronic states, such as ⁷A, ⁹A, and ¹¹A, can be formulated for the mixed-valence species [Fe^II(O₂^–)Fe^III], among which ¹¹A is found to be the ground state, while the ⁹A state lies very close to the ¹¹A state. As discussed above, in intermediate [Fe^III(O₂^2–)Fe^III], which is most likely the compound P, the O₂ molecule is coordinated simultaneously with two metal centers. The resultant species is [Fe^III(O₂^2–)Fe^III].

Our results [29] indicate that ¹¹A is the ground state; ⁹A and ⁷A lie, respectively, 6.7 and 12.0 kcal/mol above ¹¹A. Spin-density analysis shows that intermediate III in its ¹¹A state is clearly an Fe^III–Fe^III species, with approximately five spins on each Fe-center (4.87 Fe¹ and 4.89 Fe²), and some spin delocalized into the two oxygen atoms (0.19 O¹ and 0.05 O²).

Mechanism of the reaction of compound Q with methane
In this section we discuss the reaction

{Fe2(µ-O)2}Q + CH4 {Fe2(µ-O)(µ-HOCH3)}.

(3)

For this purpose we chose the smallest model of compound Q, cis-(H₂O)(NH₂)Fe(µ-O)₂(²-HCOO)₂Fe(NH₂)(H₂O), IX. Note that while spectroscopic studies show that diferryl compound Q contains two AFSC high-spin Fe^IV atoms and is diamagnetic with the exchange coupling J constant of >60 cm^–1, for simplicity we performed spin-unrestricted open-shell single-determinant calculations for the FSC high-spin states ⁹A and ¹¹A. The geometries of the reactants, intermediates, and transition states are shown in Figure 4.

Figure 4

Reactant complex, IX
As seen from Figure 4, the calculated geometric parameters and the general structural character of compound IX, which can be formally written as L₄Fe(µ-O)₂FeL₄, are consistent with experimental findings [1, 4]. Namely, it has two “µ-oxo” (or so-called “diamond”-core oxygen) atoms and two bidentate carboxylate ligands coordinated with the Fe centers. Furthermore, the diamond-core Fe₂O₂ of structures IX, as in the experimentally reported compound Q, has an asymmetric structure; one of the diamond-core O atoms is located closer to one Fe center, and the other closer to the other Fe center. The calculated diamond Fe–O distances in ⁹A are a better match to the experimental values than ¹¹A.

Methane-Q complex, X
The next step of the reaction is expected to be coordination of the incoming substrate with the compound IX. As seen in Figure 4, in general, two distinct pathways, O-side and N-side, corresponding to the substrate coordination with the bridging oxygen atoms, O¹ and O², located on the H₂O- and NH₂-sides, respectively, are possible. Similarly, according to experimental data [1], the only valid pathway is coordination of the substrate from the O-side, because of the existence of the substrate coordination pocket; coordination of the substrate from the N-side is sterically hindered and is not available. In spite of that, we investigated only the N-side pathway. We are currently studying the reaction pathway on the O-side, which will be reported separately. Here, let us suggest that the O-side pathway appears qualitatively similar to the N-side pathway.

The first step of the reaction is expected to be coordination of the incoming substrate with compound IX. Our extensive search shows that coordination of CH₄ with the diamond O-atom is more favorable than any other potentially possible coordination mode. As seen in Figure 4, coordination of the CH₄ with compound IX leads to the methane-Q complex, structure X. Since the interaction between CH₄ and structure IX is extremely weak (0.7 and 0.3 kcal/mol for the ⁹A and ¹¹A states, respectively), it is probable that complex X does not exist in reality.

The activation of the methane C–H bond is found to take place on the diamond oxygen O². Our extensive attempts to find other possible pathways for methane C–H activation either converged to the transition states XI or led to the transition states higher in energy than structures XI. As seen from Figure 4, in transition state XI the C–H bond to be broken is elongated from 1.104 Å in X to 1.271 Å and 1.296 Å for the ⁹A and ¹¹A states, respectively. Furthermore, the O–H bond is nearly formed, with 1.250 Å and 1.241 Å at the transition states, compared to 0.983 Å and 0.978 Å in the products XII for the ⁹A and ¹¹A states, respectively. These geometrical parameters clearly indicate that XI is the transition state corresponding to the H-abstraction process. The H-abstraction barriers are calculated to be 19.5 and 14.4 kcal/mol relative to the corresponding CH₄ complex X for the ⁹A and ¹¹A states, respectively. These values of the barrier are in reasonable agreement with available experimental estimates, 14–18 kcal/mol [27].

The spin densities in Table 1 for compounds XI and XII are nearly identical for ⁹A and ¹¹A, except for those on the O²..H..CH₃ fragment. For the CH₃ group, the spin densities for ⁹A and ¹¹A are –0.46 and +0.52, respectively.

Overcoming the transition state XI leads to the product XII, in which the Fe centers are bridged by one O and one OH ligand. XII is an oxo-hydroxyl complex, formally written as L₄Fe(µ-O)(µ-OH)FeL₄, with the methyl radical only weakly interacting via a C..HO interaction.

As seen in Table 1, in the intermediate XII, the CH₃ group is now a bound radical with spin densities of –0.98 and 1.00 for the ⁹A and ¹¹A states, respectively. The spin densities on Fe¹ and Fe² in XII can qualitatively be considered to correspond to Fe^IV with four spins and Fe^III with five spins, for the ⁹A and ¹¹A states, respectively.

The spin densities of the complex IX, transition states XI, and products XII presented in Table 1 can be interpreted in the following way. In structure XI, the total of eight and ten unpaired electrons (spins) of the ⁹A and ¹¹A states, respectively, are equally localized on two L_nFeO groups. According to the discussions presented above, one may consider IX (and X) in the ⁹A state as an [Fe^IV–Fe^IV] complex, and IX (and X) in the ¹¹A state as an [Fe^IV–Fe^III] mixed-valence complex. In both transition states XI, a radical center begins to develop on the CH₃ group, with spin densities of –0.46 and +0.52 for the ⁹A and ¹¹A states, respectively. In the intermediate XII, the CH₃ group becomes a bound radical with spin densities of –0.98 and 1.00 for the ⁹A and ¹¹A states, respectively. Now, the spin densities on Fe¹ and Fe² in XII can qualitatively be considered to correspond to Fe^IV and Fe^III for the ⁹A and ¹¹A states, respectively. In going from X(⁹A) to XI(⁹A), the formal oxidation state of Fe² changes from Fe^IV to Fe^III; going from X(¹¹A) (which is already Fe^III) to XI(¹¹A), no such change is required. Since the two Fe centers are coupled ferromagnetically in both the ⁹A and ¹¹A states, the spin of the CH₃ radical in both XI and XII must couple antiferromagnetically (with negative spin) and ferromagnetically (with positive spin) to make the total spin, 2S + 1, equal to 9 and 11, respectively. In the bound radical complexes XII(⁹A) and XII(¹¹A), in which the interaction of the CH₃ radical with the two iron atoms is very weak, the total energies are nearly identical, and the optimized geometries in Figure 5 are also nearly identical. In the transition states XI(⁹A) and XI(¹¹A), in which the spin is about 50% developed, the energies and geometries are similar. The fact that the system tends to adopt the same structures, whether the overall spin state 2S + 1 = 9 or 11, suggests that this mixed-spin [Fe^IV–Fe^III] electronic configuration is the preferred state of intermediate XII. The present spin-density analysis clearly demonstrates that the methane oxidation proceeds via a bound-radical mechanism, which is in good qualitative agreement with experiments of Valentine et al. [41].

Figure 5

The entire process IX + CH₄ XII is calculated to be endothermic by 11.4 and 3.0 kcal/mol for the ⁹A and ¹¹A states, respectively.

The next step in the reaction path should be the reaction of the methyl radical with the complex, which is found to occur via TS XIII. Figure 4 reveals that this transition state involves a torsion motion of the hydroxyl OH ligand before the methyl radical can add to the bridging hydroxyl ligand to form methanol. As the methyl radical approaches the hydroxyl ligand O²H from below the core Fe¹O¹Fe²O² plane, the H atom of the hydroxyl ligand leaves the coplanar tricoordinate O environment and bends upward to create a tetrahedral tetracoordinate O environment. Apparently this bend costs some energy before it is recovered by reaction of the methyl group to form a C–O covalent bond. The barrier heights for the CH₃ addition to the hydroxyl ligand calculated relative to the intermediate XII are 6.8 and 5.8 kcal/mol for the ⁹A and ¹¹A states, respectively. Obviously, this step of the reaction is not rate-determining, and can occur rather rapidly; the lifetime of the radical intermediate is expected to be too short to be detected experimentally.

Overcoming the barrier at XIII leads to the methanol complexes XIV. The overall reaction IX + CH₄ XIV is calculated to be exothermic by 31.4 and 48.6 kcal/mol for the ⁹A and ¹¹A states, respectively.

Table 1 shows that in the ¹¹A state, upon going from the XII to the XIV, Fe¹ changes its formal oxidation state from Fe^IV with four spins to Fe^III with five spins, while the spin density on the methyl radical is completely annihilated upon forming a covalent bond between CH₃ and OH. The transition state XIII has a spin distribution between that of XII and XIV. On the other hand, in the ⁹A state, upon going from the hydroxyl complex XII to the methanol complex XIV, the spin density on Fe¹ is reduced by about 0.5, corresponding to the disappearance of roughly one unpaired electron. Since Fe^V is not a stable species, it is probable that Fe¹ changes its formal oxidation state from Fe^IV with four spins to Fe^III with five formal d electrons. Because of the 2S + 1 = 9 restriction (i.e., the total number of unpaired electrons must be 8 within the Fe^III–Fe^III core), Fe¹ in XIV chose to form one d lone pair with only three spins remaining. This complex XIV(⁹A) is thus higher in energy than the corresponding complex XIV(¹¹A), in apparent violation of Hund's rule. The final step of the reaction is an elimination of the methanol molecule and regeneration of the enzyme. This step of the reaction could be a complex process; it is currently under study.

Finally, let us discuss the entire potential energy profiles of the reaction: IX + CH₄ XIV for both ⁹A and ¹¹A states, as shown in Figure 5. First, the potential energy profile does not differ much between the ⁹A and ¹¹A states. However, there are some differences in the energetics between the two states. At reactant IX, the ⁹A state is lower than the ¹¹A state by 8.8 kcal/mol, while at the first transition state XI, the energy gap between the ⁹A and ¹¹A states is reduced to 3.7 kcal/mol. The ⁹A state of IX has, from a qualitative point of view, a Fe^IV–Fe^IV core, as suggested experimentally for Q, while the ¹¹A state of IX has a Fe^IV–Fe^III core and is less stable. Once the system reaches the hydroxyl complex XII, calculations for both ⁹A and ¹¹A converge to the same electronic state with the same structure and energy, corresponding to the Fe^IV–Fe^III mixed-valence state, which interacts very weakly with the methyl radical either ferro- or antiferromagnetically. In the product methanol complex XIV, the ¹¹A state is 8.4 kcal/mol lower than the ⁹A state. Here, the preferred iron core is Fe^III–Fe^III, and each Fe has five spins, which naturally gives the ¹¹A state when ferromagnetically coupled.

Thus, the results presented above show that the MMO-catalyzed methane hydroxylation reaction proceeds via a bound-radical mechanism through the rate-determining H abstraction transition state. During the process, the oxidation state of the Fe core changes from Fe^IV–Fe^IV in IX to a mixed-valence Fe^IV–Fe^III in the short-lived intermediate XIV, and finally to Fe^III–Fe^III in XIV:IX, L₄Fe^IV(µ-O)₂Fe^IVL₄ + CH₄ XI, TS1 for CH activation (rate-determining) XII, L₄Fe^IV(µ-O)(µ-OH(· · ·CH₃))Fe^IIIL₄ (bound-radical intermediate) XIII, TS2 for CH₃ addition XIV, L₄Fe^III(OHCH₃)(µ-O)Fe^IIIL₄.

Perspective
The above calculations using small and medium-sized models have already provided new insights into the mechanisms of reactions and transformations in the MMO and R2 systems. This is because these processes are essentially electronic in nature, and the understanding of the behavior of electrons (how the spin state changes, how the substrate structure changes, and how the potential energy varies) during the processes is essential for understanding of the mechanism. However, we recognize that these simple models do not describe real metalloenzyme systems, including their entire protein environment. For instance, the existence of different amino acids on the active sites of MMOH and RNR has not been considered. To study the effects of protein environment on the catalytic activities of enzymatic systems, one must consider all (or many) of the ligands and water molecules in the real enzymes. There exists virtually no modeling calculation of enzymatic reactions in which the effects of protein environment are fully taken into account. How we accomplish this is the most serious problem to be solved, and we return to it later in Section 4.

3. The ONIOM method—A new computational tool for the study of large molecular systems

The electronic structure theory used in the above calculations for models of metalloenzyme systems is a density functional theory (DFT) method with the B3LYP functional. There exists a hierarchy of theoretical methods for calculations of the molecular structure and energy of molecular systems, as shown in Table 2.

For instance, the CCSD(T) method (coupled cluster with single and double excitations with perturbation correction for triples), one of the most reliable ab initio molecular orbital (MO) methods, can calculate the energetics, such as the enthalpy of formation, with an average error from the experiment of 2 kcal/mol or so, and provides what often is called “chemical accuracy,” the quantitative accuracy needed for chemistry. However, this method is very expensive, and the largest molecule for which this method can be routinely applied is limited at present to one containing five to seven nonhydrogen atoms. Its computer time requirement increases approximately in the order of ~N⁶, where N is the number of atoms; with this high-N dependence, the maximum size of an applicable molecule increases very slowly with the increased speed of available computers. This quantitative method cannot be applied directly to enzyme reactions in the foreseeable future. In addition, some semiquantitative methods, such as MP2/DZP or DFT(B3LYP)/DZP, give semiquantitative results with which the mechanisms of chemical reactions can be elucidated quite well, although the calculation of energies at this level of approximation is not in general good enough for determination of the absolute value of the rate constant. Some semiempirical MO methods, such as AM1 and PM3, are fast but provide semiqualitative results (i.e., conclusions derived from these methods are usually acceptable, with some exceptions). At the end of the hierarchy are the empirical classical molecular mechanics (MM) or force-field methods, such as MM3, Charmm, and Amber. These MM methods in general express molecular energy as a sum of simple functions of bond distances, bond angles, and torsion angles. With parameters fitted to experimental values, calculations are extremely rapid and give reasonable molecular and intermolecular structures and energies for systems for which the parameters are calibrated. One of the most serious drawbacks of the MM methods is that they usually cannot describe the breaking or formation of bonds (i.e., chemical reactions). Therefore, if one is interested in chemical reactions, one must use DFT (with good density functional) or MP2 level for qualitative results, and better methods for more quantitative results.

In recent years, computational chemists have become interested in chemical reactions of large molecular systems. As discussed in Section 3, for instance, we are interested in enzymatic reactions, chemical reactions in the protein environment. In other cases, one is interested in a series of organometallic reactions in homogeneous catalysis, in which the catalyst is typically a transition-metal complex with complicated ligands. In yet other cases, one is interested in chemical reactions in solution, in which some solvent molecules participate directly in the reaction or form site-specific complexes, while many other solvent molecules also affect the reaction mechanism and rate.

Using present and foreseeable future technology, it is not likely that a reliable high-level MO method can be applied to the entire molecular systems considered above. Theoreticians usually have two options. One is to adopt a simple model system to which one can apply a reasonable theoretical method. In Section 2, we did exactly this, by adopting models consisting only of metal atoms, some simplified ligands, and substrate. However, by using simplified models, one misses the effects of the protein environment, the effects of real complex ligands, or the effects of solvent molecules. These effects might change the transition-state structure and energetics, and consequently the mechanism of reaction. The second option is to adopt an inexpensive and not very reliable method with which to perform the calculation for the real (or more realistic) system. If cheap methods are used, one's results and conclusion may not be correct.

However, one can easily recognize that the chemical reaction and many other chemical problems are often a local phenomenon. A chemical reaction typically involves a few to several atoms in the action center (the core of an onion), with bond formation and breaking taking place among them. Obviously, one must use an expensive, reliable method to describe this active region of the system. The atoms in the vicinity of the action center (the second layer of the onion) participate in the reaction indirectly by pulling or pushing electrons; this region probably has to be handled with an appropriate MO method in order to take into account its electronic effects. Atoms farther away from the action center (the third layer or outer skin of the onion) are likely to provide some field, such as the long-range electrostatic interaction or nonbonding steric interaction; they can probably be treated with an inexpensive MO method or a molecular mechanics method. This kind of consideration leads to the idea of a hybrid method, a theoretical method that uses different methods for different parts of the molecular system.

The hybrid method has in fact been in use for some time. In particular, the QM/MM methods, which combine an MO method and an MM method, have been used successfully for the modeling of biological systems. These methods are usually combination-specific; i.e., a specific MO method (primarily semiempirical) and a specific MM method are combined, and additional parameters are often introduced for better agreement with experiment. We refer the readers to a recent review paper for a discussion of these methods [42].

We have recently developed the ONIOM method. At first, we used the integrated molecular orbital + molecular mechanics method (IMOMM), which integrates MO methods with the MM3 force field [43] and has a methodology slightly different from QM/MM. Soon afterward, it was realized that the extrapolation scheme in IMOMM could be generalized to combine two MO methods as well. This resulted in a combined MO + MO method, which was referred to as the integrated molecular orbital + molecular orbital method (IMOMO) [44]. Later, the integration of more than two methods was accomplished, and the entire suite of integrated methods was named the ONIOM method (our own n-layered integrated molecular orbital + molecular mechanics method) [45]. Thus, IMOMO encompasses both two-layered ONIOM2(MO:MO) and three-layered ONIOM3(MO: MO:MO), and IMOMM is in principle equivalent to ONIOM2(MO:MM) and ONIOM3(MO:MO:MM). The latest incarnation of ONIOM is available in Gaussian98, and can integrate up to three layers of most electronic structure methods and/or molecular mechanics methods from the package [46]. Various aspects of the ONIOM method have been reviewed previously [47–49]. In the present review, we briefly discuss theory and some applications of the ONIOM method.

Foundation of the ONIOM method
In the two-layered ONIOM method, the total energy of the system is obtained from three independent calculations:

EONIOM2 = Emodelhigh + Ereallow – Emodellow,

(4)

where real denotes the full system, which is treated at the low level, while model denotes the part of the system for which the energy is calculated at both high and low levels. The concept of the ONIOM method is represented schematically in Figure 6. One can see that the method can be regarded as an extrapolation scheme. Beginning at E_model^low, the extrapolation to the high-level calculation (E_model^high – E_model^low) and the extrapolation to the real system (E_real^low – E_model^low) are assumed to produce an estimate for E_real^high. For a three-layer ONIOM scheme, the energy expression can be written as

Figure 6

EONIOM3 =	Esmall modelhigh + Eintermediate modelmedium – Esmall modelmedium
	+ Ereallow–Eintermediate modellow.	(5)

The definition of the model system is rather straightforward when there is no covalent bond between the layers. If one considers a solute molecule solvated with one solvent molecule to be the real system, the solute will be the model system, and the solute–solvent interaction is naturally included in the low-level calculation of the real system. The situation becomes more complicated when a covalent bond exists between the layers. The most straightforward way to ensure that the model system is representative of the real system is to saturate the dangling bonds with link atoms, which is the scheme chosen for the ONIOM method. For example, if a methyl group is treated at the low level, a hydrogen atom can be substituted for this group in the model system. One chooses link atoms that best mimic the substituents that exist only in the real system. In practice, hydrogen atoms are good link atoms when carbon–carbon bonds are broken, but in some cases other atoms may be better [50]. One may also shift the one-electron energy of the link atom with a shift operator to mimic the electronegativity of the replaced atom [51], although we have not explored its usefulness in ONIOM in detail.

As can be seen, there is no restriction on the methods used at various levels (high, medium, low), and various MO and MM combinations discussed above can be derived. Combining different levels of MO methods is a unique feature of the ONIOM method. It turns out that MO:MO integration is rather straightforward, and virtually no special attention is required. On the other hand, the integration in ONIOM of MO and MM methods, combining two methods with very different philosophies, leads to many serious problems, as with all of the QM/MM methods. Although we do not discuss these technical problems here, we revisit this problem briefly in Section 4. The integration of two MO levels with one MM level, ONIOM3(MO:MO:MM), is unique, a feature absent from other QM/MM methods. In our MO–MM combinations, the interaction between the MO and MM regions can be treated at the MM level, i.e., with so-called mechanical embedding [52], or, alternatively, in the QM Hamiltonian, with so-called electronic embedding (in the near-future release of Gaussian software).

In the construction of the ONIOM model system, atoms that belong to the high-level layer have the same coordinates as the corresponding atoms in the real system. Even during geometry optimizations, these coordinates remain identical to one another. When no bond exists between the two layers, the first derivative of the energy with respect to the geometry is easy to obtain:

EONIOM	=	Emodelhigh	+	Ereallow	–	Emodellow	.	(6)
q		q		q		q

However, the link atoms used in the model system do not exist in the real system, and one of the main issues in this type of hybrid method is their geometrical placement, or how the geometry of the model system is related to that of the real system. We assume that the link atoms are connected to the high-level layer with the same angular and dihedral values as the link atom hosts (LAHs, the atoms replaced by the link atoms in the model system) in the real system. Now steric effects of the substituents are also taken into account in the two model system calculations. In our earlier implementation in the IMOMM method, we used fixed (standard) bond lengths between the link atoms and the high-level layer, as well as fixed bond lengths between the LAH atoms and the high-level layer. Although this scheme works well for geometry optimization, one degree of freedom is lost for each link between the high- and low-level layers, which causes problems, for example, with dynamics or frequency calculations.

In the later implementation including Gaussian98 [50], the angles and dihedrals are treated in the same manner as in the IMOMM scheme, while the bond distances between the high-level layer and the link atoms are obtained by scaling the corresponding distances between the high-level layer and the LAH atoms:

Rlink = Rhigh-level atom + g(RLAH – Rhigh-level atom),

(7)

where R_{high-level atom} denotes the atom in the high-level layer to which the link atom is connected. The scaling factor g is chosen so that reasonable bond lengths between the LAH atoms and high-level-layer atoms also yield reasonable bond lengths between the link atoms and high-level-layer atoms. In this version, the correct number of degrees of freedom ensures that the potential energy surface is well defined, and gradients and higher derivatives are available [50]. In this case, one must use the Jacobian J to convert the coordinate system for the model system to the coordinate system for the real system:

EONIOM	=	Emodelhigh	• J +	Ereallow	–	Emodellow	• J.	(8)
q		q		q		q

The Hessian or higher-order derivatives can be uniquely defined in a similar fashion. Any method for the investigation of potential energy surfaces based on conventional techniques can now be used with the ONIOM method.

Bond dissociation energies
ONIOM methods can be applied to virtually any chemical problem, but one of the first areas that has been covered in great detail is the calculation of bond dissociation energies (BDEs). In addition to several papers in which BDEs were used to illustrate the method [50, 53], a number of studies have focused exclusively on this topic. The dissociation energies of sets of larger molecules were calculated by Froese [54, 55] and Vreven¹ [56]. In both studies, a number of high and low levels were investigated, but the geometries and zero-point vibrational energies or temperature corrections were obtained at the nonintegrated B3LYP level. In most cases, a minimal model was employed (the non-H atoms only in the dissociation fragment), and the results show that the G2MS:RMP2 combination yields BDEs very close to the experimental values [G2MS [57] is our own inexpensive G2-like method, G2MS(R/U) = CCSD(T)/6-31G(d) + (R/U)OMP2/6-311 + G(2df,2p) – (R/U)OMP2/6-31G(d), which we have used extensively in our IMOMO studies]. In Table 3 we show the molecules from these studies that contain more than six non-H atoms. The error with the experimental values is usually not more than a few kcal/mol, and could probably be reduced further by fine-tuning the choice of basis sets. The primary conclusion is that a very high-level correlation method is required for the reaction center, but a low-level correlation method is still needed for the remainder of the system. Spin-contamination in the low-level method is an important problem for molecules containing aromatic substituents, and unrestricted methods have therefore proven inappropriate. We were able to predict the BDE for Ph₃C–H and Ph₃C–CH₃, for which the experimental values are not available.

Table 3	Bond dissociation energies (kcal/mol) using a minimal model. See references for specific procedures and origin of experimental data.

Reactiona Ref.b Experimental G2MS(R):RMP2c G2MS(U):RMP2c   Ph–F I 125.7 ± 2 124.6 125.2 Ph–CH3 I 101.8 ± 2 100.3 101.9 PhSiH2–H I 88.2 86.3 86.3 PhCH2–SCH3 I 61.4 ± 2 61.6 62.4 F5SO–OSF5 I 37.2 ± 2 37.3 37.7 PhCH2–H II 88.0 ± 1 90.1 90.4 MePhCH–H II 85.4 ± 1.5 87.8 88.2 Me2PhC–H II 84.4 ± 1.5 87.2 87.6 Ph2CH–H II 84 ± 2 82.6 82.9 MePh2C–H II 81 ± 2 82.6 82.9 Ph3C–H II n.a. 75.9 PhCH2–CH3 II 75.8 ± 1 79.0 79.7 MePhCH–CH3 II 74.6 ± 1.5 78.1 78.9 Me2PhC–CH3 II 73.7 ± 1.5 77.2 77.9 Ph2CH–CH3 II 72 ± 2 72.7 73.5 MePh2C–CH3 II 69 ± 2 72.9 73.6 Ph3C–CH3 II n.a. 64.1 G2MS(R):RMP2:RHFd C60:S0 T1 ( bond) III 36.1 35.1 G2MS(R):RMP2:B3LYPe Ph3C–CPh3 IV n.a. 16.4

The last two entries in Table 3, C₆₀ fullerene and hexaphenylethane (HPE), are systems too large for the G2MS:MP2 calculation, and it was necessary to introduce a third layer treated at a less expensive RHF level. In the case of C₆₀, RHF for the lowest level proved useful, and B3LYP was a good choice for HPE. The HPE example is particularly interesting; this compound is thought to be too unstable to exist [58], but our BDE value suggests that synthesis could be possible.

Diels–Alder reactions
In addition to bond dissociation energies, Diels–Alder reactions have been studied extensively with the ONIOM method. In the paper that first defined the ONIOM method, we discussed the addition of acrolein to substituted 1,3-butadienes [45], where ethylene–butadiene was used as the model system. It was expected that the electronic effect of the COH group in acrolein must also be treated with an MO method (MM is incapable of treating the electronic effects), and one was therefore included at an intermediate layer. The alkyl substituents on butadiene were not expected to play an important role, and were treated at a third, low-accuracy level. Both the optimization and energy calculations were performed with a variety of combinations. The ONIOM(B3LYP:HF:MM3) combination (the 6-31G basis set has been used throughout) yields a deviation in the activation barrier of only 1.4 kcal/mol from the B3LYP target for both the reaction of acrolein with isoprene (2-methyl-1,3-butadiene) and the reaction of acrolein with 2-tert-butyl-1,3-butadiene. The third reaction studied, ethylene + s-trans trans,trans-1,4-di-tert-butyl-1,3-butadiene, with MP4 as the target method, yields good results with the ONIOM(MP4:MP2:MM3) combination. The error in activation energy is 2.3 kcal/mol, which is significantly less than the error of 4.4 kcal/mol for MP4:HF:MM3. This is similar to the IMOMO studies of bond dissociation energies, in which we observed that results are often inadequate when the HF level is employed for a layer too close to the reaction center.

The second investigation of Diels–Alder reactions was reported in the paper that introduced G2MS [57]. The IMOMO(G2MS:MP2)//B3LYP + ZPE(B3LYP) combination was very successful, and two subsequent IMOMO studies of Diels–Alder reactions focused on this combination [59, 60]. Because G2MS calculations are too expensive for systems with more than eight non-H atoms, the method was tested by comparing the results directly to experimental values. In Table 4 we summarize the IMOMO(G2MS:MP2) results for Diels–Alder reactions for which experimental activation energies are available. The computational data is in fairly good agreement with the experimental values. In Table 5 we show the branching ratio of several reactions. The results do not always agree quantitatively with the experimental values because the branching ratio is very sensitive to small errors in our calculated energies. However, IMOMO(G2MS:MP2) reproduces qualitatively well the reversal of syn/anti reactivity for ethylene and acetylene with 1, which is caused by the repulsion between the lone pair of oxygens in 1 with the -orbitals of acetylene in the syn transition state.

Table 5	IMOMO(G2MS:MP2) and experimental (in parentheses) branching ratios in Diels–Alder reactions.

Reactants Product ratios 1,2,cis 1,2,trans 1,3,cis 1,3,trans   Acrylic acid + 2,4- pentadienoic acid 75.5(61) 22.9(22) 0.9(9) 0.7(8) syn anti Ethylene + 1 99.9(93) 0.1(7) Acetylene + 1 7.2(8) 92.8(92)

The S-value test
The most important questions when the ONIOM scheme is used are how to select the methods to be combined, and how to partition the system into high- and low-level layers. These two factors are closely related. When using the ONIOM method, one must first select the high-level method. The high-level calculation for the real system is then the target calculation that one attempts to approach with the ONIOM method. The ONIOM results converge to the target results when the low-level method approaches the high level, or when the size of the model system approaches that of the real system. Thus, there are two ways to improve the ONIOM results, both of which should be considered.

Is there any principle for selection of the appropriate computational levels and size of the model system? It depends considerably on the properties in which one is interested and on the magnitude of error one can tolerate in the ONIOM extrapolation. Unfortunately, the studies reported so far cover only a small number of chemical problems, so it is unlikely that the problem of interest to a specific reader would be among them. The first rule of thumb is that if one knows by intuition the major players and the minor players, the major players should be included in the model system to be treated at a high level, while the minor players should be treated at only a low level in the real system. For instance, if one is interested in chemical reactions in a solvated cluster (the real system), the reactants (and solvent molecules directly involved in the reaction) will constitute the model system. However, in some cases, one cannot identify the major players. For example, for the equilibrium structure of HCl(H₂O)₄ modeled by Re et al.² (Figure 7), the relative energies of one neutral complex structure (a) and two ion-pair structures (b) and (c) at ONIOM[B3LYP/D95++(d,p):BLYP/D95+(d)] (HCl at high level and water molecules at low level) do not agree with either high- or low-level ordinary MO results. A detailed analysis² shows that the interaction energy in these clusters consists of the sum or difference of large two-, three-, and multibody terms, and relatively small errors introduced by ONIOM in individual terms accumulate and easily wipe out the small (~4 kcal/mol) difference in the total binding energy among different structures.

Figure 7

We have abundant experience with the study of bond dissociation, for which it is clear which ONIOM schemes to use³ [54–56]. When an appropriate combination is not available from previous studies, chemical intuition, combined with trial and error, can be used to obtain schemes that agree with experimental data. A good example of such a procedure is given for the bond energies in C₆₀ [55]. A problem with this strategy is that the low-level method in ONIOM calculations often plays a role that is very different than if it were a stand-alone method, and it can therefore be difficult to use intuition to find the appropriate method. Furthermore, the number of possibilities (based on computational levels and partitioning) increases exponentially when the system becomes larger and an integrated method involving three or more layers is employed. A more direct way of finding an appropriate ONIOM method is preferred.

A systematic way of finding a correct ONIOM combination is by virtue of the substituent-value (S-value) test, which we have used in several of our recent ONIOM studies. Take, for instance, a relative energy E, such as bond energy or activation energy. The S-value for a certain level is defined as the difference in E between the real and the model system:

Slevel = Ereallevel