Generalized hybrid-orbital method for combining density functional theory with molecular mechanicals.

The generalized hybrid orbital (GHO) method has previously been formulated for combining molecular mechanics with various levels of quantum mechanics, in particular semiempirical neglect of diatomic differential overlap theory, ab initio Hartree-Fock theory, and self-consistent charge density functional tight-binding theory. To include electron-correlation effects accurately and efficiently in GHO calculations, we extend the GHO method to density functional theory in the generalized-gradient approximation and hybrid density functional theory (denoted by GHO-DFT and GHO-HDFT, respectively) using Gaussian-type orbitals as basis functions. In the proposed GHO-(H)DFT formalism, charge densities in auxiliary hybrid orbitals are included to calculate the total electron density. The orthonormality constraints involving the auxiliary Kohn-Sham orbitals are satisfied by carrying out the hybridization in terms of a set of L?wdin symmetrically orthogonalized atomic basis functions. Analytical gradients are formulated for GHO-(H)DFT by incorporating additional forces associated with GHO basis transformations. Scaling parameters are introduced for some of the one-electron integrals and are optimized to obtain the correct charges and geometry near the QM/MM boundary region. The GHO-(H)DFT method based on the generalized gradient approach (GGA) (BLYP and mPWPW91) and HDFT methods (B3 LYP, mPW1PW91, and MPW1 K) is tested-for geometries and atomic charges-against a set of small molecules. The following quantities are tested: 1) the C--C stretch potential in ethane, 2) the torsional barrier for internal rotation around the central C--C bond in n-butane, 3) proton affinities for a set of alcohols, amines, thiols, and acids, 4) the conformational energies of alanine dipeptide, and 5) the barrier height of the hydrogen-atom transfer between n-C4H10 and n-C4H9, where the reaction center is described at the MPW1 K/6-31G(d) level of theory.     

Testing electronic structure methods for describing intermolecular H...H interactions in supramolecular chemistry

In this article a wide variety of computational approaches (molecular mechanics force fields, semiempirical formalisms, and hybrid methods, namely ONIOM calculations) have been used to calculate the energy and geometry of the supramolecular system 2-(2'-hydroxyphenyl)-4-methyloxazole (HPMO) encapsulated in beta-cyclodextrin (beta-CD). The main objective of the present study has been to examine the performance of these computational methods when describing the short range H. H intermolecular interactions between guest (HPMO) and host (beta-CD) molecules. The analyzed molecular mechanics methods do not provide unphysical short H...H contacts, but it is obvious that their applicability to the study of supramolecular systems is rather limited. For the semiempirical methods, MNDO is found to generate more reliable geometries than AM1, PM3 and the two recently developed schemes PDDG/MNDO and PDDG/PM3. MNDO results only give one slightly short H...H distance, whereas the NDDO formalisms with modifications of the Core Repulsion Function (CRF) via Gaussians exhibit a large number of short to very short and unphysical H...H intermolecular distances. In contrast, the PM5 method, which is the successor to PM3, gives very promising results. Our ONIOM calculations indicate that the unphysical optimized geometries from PM3 are retained when this semiempirical method is used as the low level layer in a QM:QM formulation. On the other hand, ab initio methods involving good enough basis sets, at least for the high level layer in a hybrid ONIOM calculation, behave well, but they may be too expensive in practice for most supramolecular chemistry applications. Finally, the performance of the evaluated computational methods has also been tested by evaluating the energetic difference between the two most stable conformations of the host(beta-CD)-guest(HPMO) system.     

Valence bond theory for chemical dynamics

This essay provides a perspective on several issues in valence bond theory: the physical significance of semilocal bonding orbitals, the capability of valence bond concepts to explain systems with multireferences character, the use of valence bond theory to provide analytical representations of potential energy surfaces for chemical dynamics by the method of semiempirical valence bond potential energy surfaces (an early example of specific reaction parameters), by multiconfiguration molecular mechanics, by the combined valence bond-molecular mechanics method, and by the use of valence bond states as coupled diabatic states for describing electronically nonadiabatic processes (photochemistry). The essay includes both ab initio and semiempirical approaches.     

Multicoefficient Gaussian-3 calculation of the rate constant for the OH + CH4 reaction and its 12C/13C kinetic isotope effect with emphasis on the effects of coordinate system and torsional treatment

Rate constants and (12)C/(13)C kinetic isotope effects are calculated by direct dynamics for the OH + CH(4) --> H(2)O + CH(3) reaction. The electronic structure calculations required to generate the implicit potential energy surface were carried out by the high-level multicoefficient Gaussian-3/version-3 (MCG3) method and compared to two other multilevel methods, MC3BB and MC3MPW, and three density functional methods, M06-2X, BB1K, and MPW1K. The rate constants and (12)C/(13)C kinetic isotope effects are shown to depend strongly on the coordinate system used to calculate the frequencies as well as on the method used to account for the torsional anharmonicity of the lowest-frequency vibrational mode of the generalized transition states.     

Quantum dynamics of hydride transfer catalyzed by bimetallic electrophilic catalysis: synchronous motion of Mg(2+) and H(-) in xylose isomerase

Xylose isomerase exhibits a bridged-bimetallic active-site motif in which the substrate is bound to two metals connected by a glutamate bridge, and X-ray crystallographic studies suggest that metal movement is involved in the hydride transfer rate-controlling catalytic step. Here we report classical/quantal dynamical simulations of this step that provide new insight into the metal motion. The potential energy surface is calculated by treating xylose with semiempirical molecular orbital theory augmented by a simple valence bond potential and the rest of the system by molecular mechanics. The rate constant for the hydride-transfer step was calculated by ensemble-averaged dynamical simulations including both variational transition-state theory for determination of the statistically averaged dynamical bottleneck and optimized multidimensional tunneling calculations. The dynamics calculations include 25 317 atoms, with quantized vibrational free energy in 89 active-site degrees of freedom, and with 32 atoms moving through static secondary zone transition-state configurations in the quantum tunneling simulation. Our simulations show that the average Mg-Mg distance R increases monotonically as a function of the hydride-transfer progress variable z. The range of the average R along the reaction path is consistent with the X-ray structure, thus providing a dynamical demonstration of the postulated role of Mg in catalysis. We also predicted the primary deuterium kinetic isotope effect (KIE) for the chemical step. We calculated a KIE of 3.8 for xylose at 298 K, which is consistent with somewhat smaller experimentally observed KIEs for glucose substrate at higher temperatures. More than half of our KIE is due to tunneling; neglecting quantum effects on the reaction coordinate reduces the calculated KIE to 1.8.     

Theoretical investigation of hyperthermal reactions at the gas-liquid interface: O (3P) and squalane

Hyperthermal collisions (5 eV) of ground-state atomic oxygen [O ((3)P)] with a liquid-saturated hydrocarbon, squalane (C(30)H(62)), have been studied using QM/MM hybrid "on-the-fly" direct dynamics. The surface structure of the liquid squalane is obtained from a classical molecular dynamics simulation using the OPLS-AA force field. The MSINDO semiempirical Hamiltonian is combined with OPLS-AA for the QM/MM calculations. In order to achieve a more consistent and efficient simulation of the collisions, we implemented a dynamic partitioning of the QM and MM atoms in which atoms are assigned to QM or MM regions based on their proximity to "seed" (open-shell) atoms that determine where bond making/breaking can occur. In addition, the number of seed atoms is allowed to increase or decrease as time evolves so that multiple reactive events can be described. The results show that H abstraction is the most important process for all incident angles, with H elimination, double H abstraction, and C-C bond cleavage also being important. A number of properties of these reactive channels, as well as inelastic nonreactive scattering, are investigated, including angular and translational energy distributions, the effect of incident collision angle, variation with depth of the reactive event within the liquid, with the reaction site on the hydrocarbon, and the effect of dynamics before and after reaction (direct reaction versus trapping reaction-desorption).     

Analysis of classical and quantum paths for deprotonation of methylamine by methylamine dehydrogenase.

The hydrogen-transfer reaction catalysed by methylamine dehydrogenase (MADH) with methylamine (MA) as substrate is a good model system for studies of proton tunnelling in enzyme reactions--an area of great current interest--for which atomistic simulations will be vital. Here, we present a detailed analysis of the key deprotonation step of the MADH/MA reaction and compare the results with experimental observations. Moreover, we compare this reaction with the related aromatic amine dehydrogenase (AADH) reaction with tryptamine, recently studied by us, and identify possible causes for the differences observed in the measured kinetic isotope effects (KIEs) of the two systems. We have used combined quantum mechanics/molecular mechanics (QM/MM) techniques in molecular dynamics simulations and variational transition state theory with multidimensional tunnelling calculations averaged over an ensemble of paths. The results reveal important mechanistic complexity. We calculate activation barriers and KIEs for the two possible proton transfers identified-to either of the carboxylate oxygen atoms of the catalytic base (Asp428beta)-and analyse the contributions of quantum effects. The activation barriers and tunnelling contributions for the two possible proton transfers are similar and lead to a phenomenological activation free energy of 16.5+/-0.9 kcal mol(-1) for transfer to either oxygen (PM3-CHARMM calculations applying PM3-SRP specific reaction parameters), in good agreement with the experimental value of 14.4 kcal mol(-1). In contrast, for the AADH system, transfer to the equivalent OD1 was found to be preferred. The structures of the enzyme complexes during reaction are analysed in detail. The hydrogen bond of Thr474beta(MADH)/Thr172beta(AADH) to the catalytic carboxylate group and the nonconserved active site residue Tyr471beta(MADH)/Phe169beta(AADH) are identified as important factors in determining the preferred oxygen acceptor. The protein environment has a significant effect on the reaction energetics and hence on tunnelling contributions and KIEs. These environmental effects, and the related clearly different preferences for the two carboxylate oxygen atoms (with different KIEs) in MADH/MA and AADH/tryptamine, are possible causes of the differences observed in the KIEs between these two important enzyme reactions.     

Theoretical study of the dynamics of the H+CH4 and H+C2H6 reactions using a specific-reaction-parameter semiempirical Hamiltonian

We present a theoretical study of the reactions of hydrogen atoms with methane and ethane molecules and isotopomers. High-accuracy electronic-structure calculations have been carried out to characterize representative regions of the potential-energy surface (PES) of various reaction pathways, including H abstraction and H exchange. These ab initio calculations have been subsequently employed to derive an improved set of parameters for the modified symmetrically-orthogonalized intermediate neglect of differential overlap (MSINDO) semiempirical Hamiltonian, which are specific to the H+alkane family of reactions. The specific-reaction-parameter (SRP) Hamiltonian has then been used to perform a quasiclassical-trajectory study of both the H+CH4 and H+C2H6 reactions. The calculated values of dynamics properties of the H+CH4-->H2+CH3 reaction and isotopologues, including alkyl product speed distributions, diatomic product internal-state distributions, and cross sections, are generally in good agreement with experiment and with the results provided by the ZBB3 PES [Z. Xie et al., J. Chem. Phys. 125, 133120 (2006)]. The results of trajectories propagated with the SRP Hamiltonian for the H+C2H6-->H2+C2H5 reaction also agree with experiment. The level of agreement between the results calculated with the SRP Hamiltonian and experiment in both the H+methane and H+ethane reactions indicates that semiempirical Hamiltonians can be improved for not only a specific reaction but also a family of reactions.     

Variational optimization of effective atom centered potentials for molecular properties

In plane wave based electronic structure calculations the interaction of core and valence electrons is usually represented by atomic effective core potentials. They are constructed in such a way that the shape of the atomic valence orbitals outside a certain core radius is reproduced correctly with respect to the corresponding all-electron calculations. Here we present a method which, in conjunction with density functional perturbation theory, allows to optimize effective core potentials in order to reproduce ground-state molecular properties from arbitrarily accurate reference calculations within standard density functional calculations. We demonstrate the wide range of possible applications in theoretical chemistry of such optimized effective core potentials (OECPs) by means of two examples. We first use OECPs to tackle the link atom problem in quantum mechanics/molecular mechanics (QM/MM) schemes proposing a fully automatized procedure for the design of link OECPs, which are designed in such a way that they minimally perturb the electronic structure in the QM region. In the second application, we use OECPs in two sample molecules (water and acetic acid) such as to reproduce electronic densities and derived molecular properties of hybrid (B3LYP) quality within general gradient approximated (BLYP) density functional calculations.     

Mechanism of proteolysis in matrix metalloproteinase‐2 revealed by QM/MM modeling

The mechanism of enzymatic peptide hydrolysis in matrix metalloproteinase‐2 (MMP‐2) was studied at atomic resolution through quantum mechanics/molecular mechanics (QM/MM) simulations. An all‐atom three‐dimensional molecular model was constructed on the basis of a crystal structure from the Protein Data Bank (ID: 1QIB), and the oligopeptide Ace‐Gln‐Gly～Ile‐Ala‐Gly‐Nme was considered as the substrate. Two QM/MM software packages and several computational protocols were employed to calculate QM/MM energy profiles for a four‐step mechanism involving an initial nucleophilic attack followed by hydrogen bond rearrangement, proton transfer, and C? N bond cleavage. These QM/MM calculations consistently yield rather low overall barriers for the chemical steps, in the range of 5–10 kcal/mol, for diverse QM treatments (PBE0, B3LYP, and BB1K density functionals as well as local coupled cluster treatments) and two MM force fields (CHARMM and AMBER). It, thus, seems likely that product release is the rate‐limiting step in MMP‐2 catalysis. This is supported by an exploration of various release channels through QM/MM reaction path calculations and steered molecular dynamics simulations. © 2015 Wiley Periodicals, Inc.     

Sensitivity of the H + CH3 → CH4 recombination rate constant to the shape of the CH stretching potential

Using a potential-energy surface obtained in part from ab initio calculations, the H + CH₃ → CH₄ bimolecular rate constant at T = 300 K is determined from a Monte Carlo classical trajectory study. Representing the CH stretching potential with a standard Morse function instead ofthe ab initio curve increases the calculated rate constant by an order of magnitude. The experimental recombination rate constant is intermediate of the rate constants calculated with the Morse and ab initio stretching potentials.Two properties of the H + CH₃ α CH₄ potential-energy surface which significantly affect the recombination rate constant are the shape of the CH stretching potential and the attenuation of the H₃CH bending frequencies. Ab initio calculations with a hierarchy of basis sets and treatment of electron correlation indicate the latter is properly described [13]. The exact shape of the CH stretching potential is not delineated by the ab initio calculations, since the ab initio calculations are not converged for bond lengths of 2.0–3.0 Å [12]. However, the form of this stretching potential deduced from the highest-level ab initio calculations, and fit analytically by eq. (2), is significantly different from a Morse function. The experimental recombination rate constant is intermediate of the rate constants calculated with the Morse and ab initio CH stretching potentials. This indicates that the actual CH potential energy curve lies between the Morse and ab initio curves. This is consistent with the finding that potential energy curves for diatomics are not well described by a Morse function [12].     

Investigation of the mechanism of the cell wall DD-carboxypeptidase reaction of penicillin-binding protein 5 of Escherichia coli by quantum mechanics/molecular mechanics calculations

Penicillin-binding protein 5 (PBP 5) of Escherichia coli hydrolyzes the terminal D-Ala-D-Ala peptide bond of the stem peptides of the cell wall peptidoglycan. The mechanism of PBP 5 catalysis of amide bond hydrolysis is initial acylation of an active site serine by the peptide substrate, followed by hydrolytic deacylation of this acyl-enzyme intermediate to complete the turnover. The microscopic events of both the acylation and deacylation half-reactions have not been studied. This absence is addressed here by the use of explicit-solvent molecular dynamics simulations and ONIOM quantum mechanics/molecular mechanics (QM/MM) calculations. The potential-energy surface for the acylation reaction, based on MP2/6-31+G(d) calculations, reveals that Lys47 acts as the general base for proton abstraction from Ser44 in the serine acylation step. A discrete potential-energy minimum for the tetrahedral species is not found. The absence of such a minimum implies a conformational change in the transition state, concomitant with serine addition to the amide carbonyl, so as to enable the nitrogen atom of the scissile bond to accept the proton that is necessary for progression to the acyl-enzyme intermediate. Molecular dynamics simulations indicate that transiently protonated Lys47 is the proton donor in tetrahedral intermediate collapse to the acyl-enzyme species. Two pathways for this proton transfer are observed. One is the direct migration of a proton from Lys47. The second pathway is proton transfer via an intermediary water molecule. Although the energy barriers for the two pathways are similar, more conformers sample the latter pathway. The same water molecule that mediates the Lys47 proton transfer to the nitrogen of the departing D-Ala is well positioned, with respect to the Lys47 amine, to act as the hydrolytic water in the deacylation step. Deacylation occurs with the formation of a tetrahedral intermediate over a 24 kcal x mol(-1) barrier. This barrier is approximately 2 kcal x mol(-1) greater than the barrier (22 kcal x mol(-1)) for the formation of the tetrahedral species in acylation. The potential-energy surface for the collapse of the deacylation tetrahedral species gives a 24 kcal x mol(-1) higher energy species for the product, signifying that the complex would readily reorganize and pave the way for the expulsion of the product of the reaction from the active site and the regeneration of the catalyst. These computational data dovetail with the knowledge on the reaction from experimental approaches.     

Hydride transfer catalyzed by xylose isomerase: mechanism and quantum effects

We have applied molecular dynamics umbrella-sampling simulation and ensemble-averaged variational transition state theory with multidimensional tunneling (EA-VTST/MT) to calculate the reaction rate of xylose-to- xylulose isomerization catalyzed by xylose isomerase in the presence of two Mg2+ ions. The calculations include determination of the free energy of activation profile and ensemble averaging in the transmission coefficient. The potential energy function is approximated by a combined QM/MM/SVB method involving PM3 for the quantum mechanical (QM) subsystem, CHARMM22 and TIP3P for the molecular mechanical (MM) environment, and a simple valence bond (SVB) local function of two bond distances for the hydride transfer reaction. The simulation confirms the essential features of a mechanism postulated on the basis of kinetics and X-ray data by Whitlow et al. (Whitlow, M.; Howard, A. J.; Finzel, B. C.; Poulos, T. L.; Winborne, E.; Gilliland, G. L. Proteins 1991, 9, 153) and Ringe, Petsko, and coworkers (Labie, A.; Allen, K.-N.; Petsko, G. A.; Ringe, D. Biochemistry 1994, 33, 5469). This mechanism involves a rate-determining 1,2-hydride shift with prior and post proton transfers. Inclusion of quantum mechanical vibrational energy is important for computing the free energy of activation, and quantum mechanical tunneling effects are essential for computing kinetic isotope effects (KIEs). It is found that 85% of the reaction proceeds by tunneling and 15% by overbarrier events. The computed KIE for the ratio of hydride to deuteride transfer is in good agreement with the experimental results. The molecular dynamics simulations reveal that proton and hydride transfer reactions are assisted by breathing motions of the mobile Mg2+ ion in the active site, providing evidence for concerted motion of Mg2+ during the hydride transfer step.     

Ab initio and direct quasiclassical-trajectory study of the F+CH4-->HF+CH3 reaction

We present an electronic structure and dynamics study of the F+CH4-->HF+CH3 reaction. CCSD(T)/aug-cc-pVDZ geometry optimizations, harmonic-frequency, and energy calculations indicate that the potential-energy surface is remarkably isotropic near the transition state. In addition, while the saddle-point F-H-C angle is 180 degrees using MP2 methods, CCSD(T) geometry optimizations predict a bent transition state, with a 153 degrees F-H-C angle. We use these high-quality ab initio data to reparametrize the parameter-model 3 (PM3) semiempirical Hamiltonian so that calculations with the improved Hamiltonian and employing restricted open-shell wave functions agree with the higher accuracy data. Using this specific-reaction-parameter PM3 semiempirical Hamiltonian (SRP-PM3), we investigate the reaction dynamics by propagating quasiclassical trajectories. The results of our calculations using the SRP-PM3 Hamiltonian are compared with experiments and with the estimates of two recently reported potential-energy surfaces. The trajectory calculations using the SRP-PM3 Hamiltonian reproduce quantitatively the measured HF vibrational distributions. The calculations also agree with the experimental HF rotational distributions and capture the essential features of the excitation function. The results of the SRP semiempirical Hamiltonian developed here clearly improve over those using the two prior potential-energy surfaces and suggest that reparametrization of semiempirical Hamiltonians is a promising strategy to develop accurate potential-energy surfaces for reaction dynamics studies of polyatomic systems.     

Theoretical studies of hyperthermal O(3P) collisions with hydrocarbon self-assembled monolayers

We present a dynamics study of inelastic and reactive scattering processes in collisions of hyperthermal (5 eV) O(3P) atoms with a hydrocarbon self-assembled monolayer (SAM). Molecular-dynamics simulations are carried out using a quantum mechanics/molecular mechanics (QM/MM) interaction potential that uses a high quality semiempirical Hamiltonian for the QM part and the MM3 force field for the MM part. A variety of products coming from reaction are identified, including H abstraction to generate OH, O atom addition to the SAM with subsequent elimination of H atoms, and direct C-C breakage. The C-C breakage mechanism provides a pathway for significant surface mass loss in single reactive events whereas the O addition-H elimination channel leads to surface oxidation. Reaction probabilities, product energy, and angular distributions are examined to gain insight on polymer erosion in low Earth orbit conditions and on fundamentals of inelastic and reactive hyperthermal gas-surface interactions.     

Benchmark calculations on models of the phosphoryl transfer reaction catalyzed by protein kinase A

Protein phosphorylation has been proved to be of great importance in many stages of cell life. In the last few years, its reaction mechanism has been extensively studied. In this work we present the analysis of the performances of several computational methods with different computational costs (from multilevel to semiempirical) to point out the best method to be used at each level in the study of phosphoryl transfer. Finally, we center on the semiempirical methods, and mainly on the AM1/d Hamiltonian with different sets of parameters, which will permit hybrid quantum mechanics/molecular mechanics (QM/MM) free energy calculations on big models at an acceptable computational cost. We have used quite a large set of molecules and model reactions to test the computational methods, reproducing all the chemical steps involved in the mainly accepted reaction pathways for the protein phosphorylation. In the end, we also present the results for an enlarged model, cut out from an entire biological model: we compare the 2-D PES at the B3LYP and AM1/d levels with the purpose of obtaining a correction for the semiempirical method. The AM1/d-PhoT semiempirical parameterization corrected using single-point energy calculations at the B3LYP/MG3S level seems to be suitable to carry out reliable QM/MM calculations of the complete biological system.     

Kinetics study of the OH + alkene --> H2O + alkenyl reaction class

In this paper we report on the kinetics of hydrogen abstraction for the OH + alkene reaction class, using the reaction class transition state theory (RC-TST) combined with the linear energy relationship (LER) and the barrier height grouping (BHG) approaches. Parameters for the RC-TST were derived from theoretical calculations using a set of 15 reactions representing the hydrogen abstractions from the terminal and nonterminal carbon sites of the double bond of alkene compounds. Both the RC-TST/LER, where only reaction energy is needed at either density functional theory BH&HLYP or semiempirical AM1 levels, and RC-TST/BHG, where no additional information is required, are found to be promising methods for predicting rate constants for a large number of reactions in this reaction class. Detailed error analyses show that, when compared to explicit theoretical calculations, the averaged systematic errors in the calculated rate constants using both the RC-TST/LER and RC-TST/BHG methods are less than 25% in the temperature range 300-3000 K. The estimated rate constants using these approaches are in good agreement with available data in the literature.     

Ab initio and direct quasiclassical-trajectory study of the Cl + CH4-->HCl + CH3 reaction

We present an electronic structure and dynamics study of the Cl + CH(4)--> HCl + CH(3) reaction. We have characterized the stationary points of the ground-state potential-energy surface using various electronic structure methods and basis sets. Our best calculations, CCSD(T) extrapolated to the complete basis-set limit based on geometries and harmonic frequencies obtained at the CCSD(T)/aug-cc-pvtz level, are in agreement with the experimental reaction energy and indirect measurements of the barrier height. Using ab initio information, we have reparametrized a semiempirical Hamiltonian so that the predictions of the improved Hamiltonian agree with the higher-level calculations in various regions of the potential-energy surface. This improved semiempirical Hamiltonian is then used to propagate quasiclassical trajectories and characterize the reaction dynamics. The good agreement of the calculated HCl rotational and angular distributions with the experiment indicates that reparametrizing semiempirical Hamiltonians is a promising approach to derive accurate potential-energy surfaces for polyatomic reactions. However, excessive energy leakage from the initial vibrational energy of the CH(4) molecule to the reaction coordinate in the trajectory calculations calls into question the suitability of the standard quasiclassical-trajectory method to describe energy partitioning in polyatomic reactions.     

Quantum mechanical/molecular mechanical study on the mechanism of the enzymatic Baeyer-Villiger reaction

We report a combined quantum mechanical/molecular mechanical (QM/MM) study on the mechanism of the enzymatic Baeyer-Villiger reaction catalyzed by cyclohexanone monooxygenase (CHMO). In QM/MM geometry optimizations and reaction path calculations, density functional theory (B3LYP/TZVP) is used to describe the QM region consisting of the substrate (cyclohexanone), the isoalloxazine ring of C4a-peroxyflavin, the side chain of Arg-329, and the nicotinamide ring and the adjacent ribose of NADP(+), while the remainder of the enzyme is represented by the CHARMM force field. QM/MM molecular dynamics simulations and free energy calculations at the semiempirical OM3/CHARMM level employ the same QM/MM partitioning. According to the QM/MM calculations, the enzyme-reactant complex contains an anionic deprotonated C4a-peroxyflavin that is stabilized by strong hydrogen bonds with the Arg-329 residue and the NADP(+) cofactor. The CHMO-catalyzed reaction proceeds via a Criegee intermediate having pronounced anionic character. The initial addition reaction has to overcome an energy barrier of about 9 kcal/mol. The formed Criegee intermediate occupies a shallow minimum on the QM/MM potential energy surface and can undergo fragmentation to the lactone product by surmounting a second energy barrier of about 7 kcal/mol. The transition state for the latter migration step is the highest point on the QM/MM energy profile. Gas-phase reoptimizations of the QM region lead to higher barriers and confirm the crucial role of the Arg-329 residue and the NADP(+) cofactor for the catalytic efficiency of CHMO. QM/MM calculations for the CHMO-catalyzed oxidation of 4-methylcyclohexanone reproduce and rationalize the experimentally observed (S)-enantioselectivity for this substrate, which is governed by the conformational preferences of the corresponding Criegee intermediate and the subsequent transition state for the migration step.