|Home | About | Journals | Submit | Contact Us | Français|
Application of combined quantum and molecular mechanical (QM/MM) methods focuses on predicting activation barriers and the structures of stationary points for organic and enzymatic reactions. Characterization of the factors that stabilize transition structures in solution and in enzyme active sites provides a basis for design and optimization of catalysts. Continued technological advances allowed expansion from prototypical cases to mechanistic studies featuring detailed enzyme and condensed-phase environments with full integration of the QM calculations and configurational sampling. This required improved algorithms featuring fast QM methods, advances in computing changes in free energies including free-energy perturbation (FEP) calculations, and enhanced configurational sampling. In particular, the present Account highlights development of the PDDG/PM3 semiempirical QM method, computation of multidimensional potentials of mean force (PMF), incorporation of on-the-fly QM in Monte Carlo (MC) simulations, and a polynomial quadrature method for efficient modeling of proton-transfer reactions.
The utility of this QM/MM/MC/FEP methodology is illustrated for a variety of organic reactions including substitution, decarboxylation, elimination, and pericyclic reactions. Comparison with experimental kinetic results on medium effects has verified the accuracy of the QM/MM approach in the full range of solvents from hydrocarbons to water to ionic liquids. Corresponding results from ab initio and density functional theory (DFT) methods with continuum-based treatments of solvation reveal deficiencies, particularly for protic solvents. Also summarized in this Account are three specific QM/MM applications to biomolecular systems: (1) a recent study that clarified the mechanism for the reaction of 2-pyrone derivatives catalyzed by macrophomate synthase as a tandem Michael-aldol sequence rather than a Diels-Alder reaction; (2) elucidation of the mechanism of action of fatty acid amide hydrolase (FAAH), an unusual Ser-Ser-Lys proteolytic enzyme; and, (3) the construction of enzymes for Kemp elimination of 5-nitrobenzisoxazole that highlights the utility of QM/MM in the design of artificial enzymes.
QM/MM methods permit the modeling of bond-making and bond-breaking processes, and treatment of systems much larger than by QM alone and/or that lack MM parameters.1–6 Typically, a central region is described by QM methods and the remainder of the system is represented by MM. Many alternatives are possible depending on choices for the QM, MM, representation of the solvent, QM/MM interface, and configurational sampling.2 An early approach for organic reactions in solution began by obtaining a minimum-energy reaction path from ab initio calculations in the gas phase. The effects of solvation were then determined for this pathway with importance sampling or free-energy perturbation (FEP) calculations using Metropolis Monte Carlo (MC) simulations.6 In their day, such calculations were massive, typically requiring extensive software development and access to what were then supercomputers, like the Cray 1 and Cyber 205. Application to substitutions, additions, and pericyclic reactions provided fascinating insights on variations in solvation, especially hydrogen bonding, along the reaction paths.6 Nevertheless, clear limitations of this “QM + MM” approach are the use of the same reaction path in different media and the lack of explicit polarization effects. In time, the methodology progressed to eliminate the reaction-path constraints and to integrate on-the-fly QM calculations with the liquid-state simulations.2,7 In our approach, the solutes or key parts of the reacting system are treated by QM, the solvent is part of the MM region, solute-solvent interactions consist of Coulombic terms using the QM and MM charges and Lennard-Jones (LJ) interactions using standard MM parameters, and Monte Carlo statistical mechanics provides the configurational sampling.7a This Account summarizes the development of this QM/MM/MC methodology including the PDDG/PM3 semiempirical QM method along with recent applications to organic and enzymatic reactions.
The reacting systems are surrounded by 500–1000 solvent molecules, which are included explicitly using the OPLS-AA force field for non-aqueous solvents8 and the TIP4P model for water.9 The QM/MM computations are carried out with the BOSS and MCPRO programs for organic and enzymatic reactions, respectively.10 Periodic boundary conditions are used for the organic reactions, while the active sites of enzymes are embedded in water caps. Solute-solvent and solvent-solvent interactions are typically truncated using 12-Å residue-based cutoffs with smoothing over the last 0.5 Å. The energy and wave function for the QM region are obtained from single-point calculations upon each attempted MC move of a QM element. The non-bonded potential energy between the QM and MM regions is given by Coulomb and Lennard-Jones interactions (eq 1). The σ and ε parameters are taken from the OPLS-AA force field along with the atomic charge qi for MM atoms.
An important issue is the computation of charges for the QM atoms. Our preference has been to use the Cramer-Truhlar CMx models, which have been optimized for reproduction of gas-phase dipole moments.11 Based on tests for performance in solution, our choice for neutral molecules has been CM1A (AM1-based) or CM3P (PM3-based) charges scaled by ca. 1.14.12 The 1.14*CM1A charges yield a mean unsigned error (mue) of 1.0 kcal/mol for free energies of hydration, while the mue is 0.7 kcal/mol with OPLS-AA charges.12 These choices provide a balanced treatment such that acceptable results are consistently obtained for free energies of hydration and medium effects on equilibria and reaction rates.7,13–19
For proper study of organic and enzymatic reactions, it is essential that extensive sampling of the substrate, protein, and solvent molecules be carried out to obtain configurationally averaged free-energy results. Without adequate sampling, even ab initio QM/MM methods can yield spurious findings.20 However, adequate sampling is accompanied by the need to perform large numbers of QM calculations. For example, a recent QM/MM study of Diels-Alder reactions in solution, a relatively simple case, required 3.5 million single-point QM calculations to obtain acceptable convergence for each one-dimensional (1-D) free-energy profile.21 Thus, highly efficient QM methods are needed. In our case, semiempirical QM (SQM) methods have been explored, though there are alternative approaches1–5 including the empirical valence bond (EVB) and ONIOM methods.1,4
The most common SQM choices, AM122 and PM3,23 were developed in the 1980s and followed the MNDO24 formalism. Subsequent attempts to re-optimize the parameter sets using modern optimizers provided little improvement.25,26 However, SQM methods suffer from several well-known problems.22,23 Common errors, which also plague ab initio calculations using basis sets without polarization functions, occur for molecules with small rings or adjacent heteroatoms. For the SQM methods, blame could be attributable to both the valence-only sp basis sets and neglect of differential overlap (NDO). Ameliorating efforts included reintroduction of the overlap matrix into the secular equations in MNDO yielding NO-MNDO, and the addition of d-orbitals on the first row with MNDO and AM1.26 The mean unsigned error (mue) for heats of formation was reduced from 8.4 kcal/mol using MNDO to 6.8 kcal/mol with NO-MNDO for a diverse set of 622 neutral, closed-shell molecules containing C, H, N, and O atoms, and it also gave modest improvements for the activation barriers of nine pericyclic reactions.26 The addition of d-orbitals for 217 hydrocarbons excluding alkynes, reduced the mue for ΔHf with MNDO from 9.3 to 6.8 and for AM1 from 5.3 to 3.8 with ring compounds improving from 7.1 to 4.5 kcal/mol. However, the error for alkynes was unacceptably large. The problem was traced to overly repulsive 2-center resonance integrals for short d-d interactions, which is fixable, in principle, with a damping function.
In the most successful effort, improvement of the core repulsion formula (CRF) was considered. The principal difference between AM1 and MNDO is the addition of multiple Gaussians to the MNDO CRF, e.g., 7 Gaussians for a C–H interaction in AM1.22. Starting with MNDO, we instead added the pairwise-distance-directed Gaussian expression in eq 2, which uses 4 terms for an A–B interaction and 3 for A-A. This method, PDDG/MNDO, has fewer parameters than AM1, but upon optimization reduced the ΔHf mue for the 622 molecule set from 8.4 (MNDO) and 6.7 (AM1) to 5.2 kcal/mol.25 The mue with PM3 for the 622 compounds is only 4.4; however, this was surpassed by adding the PDDG terms to PM3 yielding PDDG/PM3 and a mue of 3.2.25 Performance with the alternative SQM method, SCC-DFTB, was found to be intermediate between AM1 and PM3.27 With PDDG/PM3, well-known problems with homologation and hydrocarbon branching were corrected, and results for diverse isomerization energies were found to be more accurate (mue = 1.6) than from B3LYP/6-31G* (mue = 2.3). ΔHf results and, therefore, heats of reaction are significantly improved with PDDG/PM3 over B3LYP-based DFT methods; e.g., for the G3 dataset, the mues are 3.2 for PDDG/PM3 and 7.2 for B3LYP/6-311+G(3df,2p). Notably, the listed SQM methods all reproduce MP2/cc-pVTZ molecular geometries with average errors for bond lengths, bond angles, and dihedral angles of only ca. 0.01 Å, 1.5°, and 3°.27 However, a general weak point for SQM methods is the description of hydrogen-bonding; this issue is largely avoided in our QM/MM studies by the MM representation of the solvent and use of the CM1A charges for QM atoms.
The PDDG methods were extended to other elements.28 Dramatic improvements were obtained for 422 halogen-containing molecules with ΔHf mues for MNDO (14.0), AM1 (11.1), PM3 (8.1), PDDG/MNDO (6.6), and PDDG/PM3 (5.6).28a The error here for PDDG/PM3 comes predominantly from polyfluoro compounds. Hydrogen-bonded systems, ion-molecule complexes, and transition states were included and yielded results similar to those from B3LYP/6-311++G(d,p). Barrier heights for SN2 reactions involving Cl, Br, and I are particularly well described. PDDG/PM3 parameterization for Si, P, and S was also reported with substantial improvements for sulfur.28b The PDDG/PM3 results for hypervalent S-containing compounds are significantly better than from MNDO/d and bond-length errors were reduced by 50% over PM3. Finally, the performance of B3LYP-based DFT methods was addressed for the full set of 622 C, H, N, and O-containing molecules with optimization of atomic-energy offsets for the ΔHf calculations as with the SQM methods.29 At the highest level tested, B3LYP/6-31+G(d,p), the accuracy for heats of formation and isomerization energies is essentially identical to that with the far simpler and faster PDDG/PM3 method.29
Accurate calculation of free-energy changes is also needed for the characterization of chemical processes. FEP calculations are based on the Zwanzig equation (eq 3), which relates the free energy difference between an initial (0) and final (1) state to a configurational average involving their potential energy difference.6,30 For application of eq 3 to the medium effect on the interconversion of two molecular entities, A and B, the thermodynamic cycle in Scheme 1 is used. The A to B FEP conversion is performed in both media, and the medium effect is given by eq 4. For acceptable convergence of eq 3, the mutation is normally broken into a series of steps or “windows” that is characterized by a coupling parameter λ that ranges from 0 to 1 as A goes to B. The total ΔG is then the sum over the contributions from each window.
For reactions, free energy changes are normally computed as a function of an inter- or intra-molecular reaction coordinate, e.g., an interatomic distance or dihedral angle. The resultant free-energy profile is a “potential of mean force” (PMF). FEP calculations are performed between adjacent points on the reaction coordinate and the results are combined to yield the PMF. The geometrical change between points must be small, e.g., Δr = 0.01–0.05 Å, to ensure again acceptable convergence. It is also common to obtain multi-dimensional PMFs by perturbing between points on a grid defined by multiple reaction coordinates.5,31 For most cases covered here, 2-D PMFs were computed using the lengths of key making and breaking bonds as the reaction coordinates.
Substitution reactions have long been a testing ground for QM/MM methods.6b,32 Indeed, condensed-phase studies with the present procedures were first carried out for SN2 reactions of chloride ion with methyl, ethyl, and neopentyl chlorides (Figure 1). The transition states were located by mapping the free energy as a function of the two C-Cl distances and the Cl-C-Cl angle.33 PDDG/PM3 and CBS-QB3 results for the gas-phase reactions were virtually identical, and all aspects of the ΔG‡ predictions in solution were in good accord with experiment: the absolute ΔG‡ values, the barrier lowering in dipolar aprotic solvents, and the effects of the electrophile’s structure on reactivity. The rate retardation for SN2 reactions on neopentyl halides was confirmed to be a steric effect in all media.,33 Modeling of SN2 methyl-transfer reactions of sulfonium and ammonium salts in aqueous solution also proceeded well using a combination of CBS-QB3 and PDDG-PM3/QM/MM calculations.34
Similarly, the ΔG‡ values for the SNAr reaction of N3− + p-FPhNO2 from PDDG/PM3-based QM/MM calculations and experiment generally matched well including the rate increase by a factor of ca. 104 upon transfer from methanol to DMSO (Table 1).13 The computed activation barrier in water was relatively overestimated. This first QM/MM study of an SNAr reaction provided detailed structural characterization of the transition states and intermediate Meisenheimer complexes in four solvents. Parallel B3LYP/6-311+G(2d,p) calculations with solvation treated by the polarizable continuum model (PCM)35 did not reproduce the observed solvent effects (Table 1). The explicit description of hydrogen bonding is necessary.
For Menshutkin reactions, H3N + CH3X → H3NCH3+ + X−, PDDG/PM3 gas-phase ΔHs of 124, 117, and 102 kcal/mol compare well with the experimental values of 123, 115, and 107 kcal/mol for X = Cl, Br, and I.5 Subsequent QM/MM calculations for X = Cl in water and DMSO gave ΔG‡ values of 25.8 and 30.1 kcal/mol, which are consistent with prior computational results and an experimental value of 23.5 kcal/mol for the X = I reaction in water.5 In contrast to the SN2 and SNAr reactions with anionic nucleophiles, the Menshutkin reactions are accelerated in water in view of the charge development.32e Structural results also showed the expected earlier TS in water, r(CN) and r(CCl) of 2.14 and 2.18 Å, compared to 2.02 and 2.20 Å in DMSO, and gas-phase values of 1.82 and 2.35 Å. Interestingly, the bond length changes for this type-II SN2 reaction are similar to those for the type-I reactions, e.g., for Cl− + EtCl, r(CCl) for the TS increases from 2.29 Å in the gas phase to 2.47 Å in water. This challenges the traditional view that changes in solvent should not lead to changes in TS structure for the type-I reactions.37
QM/MM methodology has also applied to Kemp decarboxylations of 3-carboxy-benzisoxazoles,14 decarboxylation of N-carboxy-2-imidazolidinone (a biotin model),15 and Cope eliminations of threo- and erythro-N,N-dimethyl-3-phenyl-2-butylamine oxide.16 These reactions experience large rate accelerations upon transfer from protic to aprotic solvents. Both the relative and absolute ΔG‡ values were well reproduced with PDDG/PM3-based QM/MM simulations; illustrative computed (±0.3 kcal/mol) and experimental (±1.0 kcal/mol) results are listed in Table 2 for the Kemp reaction of 1 and for decarboxylation of the biotin model.14c,15b In the substituted case 3, an internal hydrogen bond is responsible for slow elimination with little solvent dependence.38 Analogous QM/MM calculations with AM1 did not reproduce this effect owing to poor representation of the hydrogen bond, while PDDG/PM3 performed well for 1–3.14c Furthermore, it was demonstrated that ion-pairing with the tetramethylguanidinium (TMG) counterion must be considered to model correctly the Kemp reaction in aprotic solvents like chloroform (Figure 2).14c The computed ΔG‡ values are 27.7 and 14.8 kcal/mol for reaction with and without TMG, while the experimental ΔG‡ is 24.0 kcal/mol.38
Ab initio and DFT calculations coupled with continuum solvation were again found to underestimate solvent effects, as illustrated by the PCM-based results in Table 3 for the Cope elimination.16 Furthermore, simulation of mixed solutions is ill-defined using continuum models. However, QM/MM/MC remains robust, as found for the biotin decarboxylation (Figure 3); e.g., computed ΔG‡s of 23.6, 22.8, and 18.5 kcal/mol were obtained in 33%, 66%, and 100% v/v acetonitrile/water mixtures compared to 22.7, 22.0, and 19.5 kcal/mol from experiment.15b There is clear de-mixing of the solvents to maximize hydrogen bonding to the solute; a little water goes a long way, so substantial rate acceleration is not obtained until the pure acetonitrile limit is approached.
Finally, for a case where there is more mechanistic uncertainty, the uncatalyzed decomposition of urea in water was examined.41 Elimination of ammonia after water-assisted rearrangement to the H3NCONH zwitterion was found to be preferred over hydrolysis by an addition/elimination mechanism. The computed barrier, 37 kcal/mol, is somewhat higher than experimental estimates of ca. 30 kcal/mol, while the 3-kcal/mol preference for the ammonia elimination is consistent with an experimental rate ratio of ca. 150 for the two processes.
Three Diels-Alder reactions, cyclopentadiene with acrylonitrile, methyl vinyl ketone, and 1,4-napthoquinone (Figure 3), were initially investigated in water using AM1/MM and by computing a 1–D PMF as a function of the distance between the midpoint of C1 and C4 for the diene and the center of the C=C bond for the dienophile.21 Though the absolute ΔG‡ values were overestimated, the effects of hydration were accurately reproduced. The same approach was followed in a QM/MM study of the reaction of cyclopentadiene and methyl acrylate in two room-temperature ionic liquids, 1-ethyl-3-methylimidazolium (EMIM) tetrachloroaluminate and heptachlorodialuminate.18 The computed ΔΔG‡ (kcal/mol) values of 0.0, 0.6 and −2.9 in water, EMIM-AlCl4, and EMIM-Al2Cl7 again mirrored well the experimental values of 0.0, 0.52, and −1.36. A recent reinvestigation of the three Diels-Alder reactions featured QM/MM/MC computation of full 2-D PMFs as a function of both forming bond lengths using PDDG/PM3 for the QM.17 The 2-D approach enables reliable location of the transition states in an automated manner. The ΔΔG‡ values were reproduced well in four solvents, while corresponding ab initio calculations with a continuum solvation model failed to yield the large rate increases in water (Table 4). Early work established the importance of enhanced hydrogen bonding at the transition states for catalysis of pericyclic reactions in protic solvents;6c the concept continues to be actively pursued in developing chiral hydrogen-bond donors as asymmetric catalysts.43
Two ene reactions, tetramethylethylene with singlet oxygen and 4-phenyl-1,2,4-triazoline-3,5-dione (PTAD), were also modeled.19,31 The PDDG/PM3-based QM/MM/FEP calculations found the mechanisms to be medium-dependent with different intermediates in the PTAD system and with conversion from a concerted to stepwise mechanism when proceeding from the gas phase to solution for the 1O2 case. Reasonable free energies of activation were obtained, e.g., 14.9 kcal/mol in acetonitrile (exptl. = 15.0 kcal/mol) for the PTAD system, and 12.5 kcal/mol in cyclohexane (exptl. = 8.8 kcal/mol) for 1O2. However, DFT/CPCM calculations fared poorly predicting a near-constant barrier of 25 kcal/mol in diverse solvents for the PTAD case.
Initial coordinates are obtained from crystal structures for protein-ligand complexes. For large proteins, 150–200 residues nearest the active site are retained. A utility program is used to convert the raw PDB file into a Z-matrix suitable for MCPRO by adding missing hydrogens, performing the residue truncation and capping, and assigning protonation states and OPLS-AA atom types.10 The substrate is placed in the active site and a ca. 25-Å radius water cap containing ca. 1000 water molecules is added. The MC simulations then involve attempted moves of protein side chains, optionally the backbone, substrates, water molecules, and any ions. The MC run for each FEP window consists of ca. 20 million configurations of equilibration and 50 million configurations of averaging. Attractive features of MC simulations for proteins include ease of implementation of geometrical constraints, ease of performance of FEP calculations without endpoint instabilities, efficient sampling in internal coordinates, and uniform, exact temperature and pressure control.10,44 Enzyme residues that are intimately involved in the reaction are included in the QM region. The interface between the QM and MM regions often involves the use of “link atoms”, which can be implemented in multiple ways.2 In our approach, the QM residues are capped with hydrogens as link atoms placed on adjacent Cα atoms.45 The QM-MM connection also requires the inclusion of classical bond-stretching, angle-bending, and torsion terms for atoms at the interface.
Macrophomate synthase (MPS) catalyzes the transformation of 2-pyrone derivatives to benzoate analogues in fungi. The transformation involves three reactions: decarboxylation of oxalacetate to produce pyruvate enolate, two C-C bond formations between a 2-pyrone and pyruvate enolate that form a bicyclic intermediate, and final decarboxylation with concomitant dehydration. Although the second step might be a long-sought enzymatic Diels-Alder reaction, proof that the C-C bond formations are concerted is lacking. A tandem Michael-aldol sequence is an alternative (Scheme 2). Report of a crystal structure for MPS complexed with pyruvate enolate and Mg2+ (see Conspectus Figure)46 enabled QM/MM/MC/FEP investigation of the two pathways.45 1–D PMFs were computed for the Michael and aldol steps, which readily located the illustrated Michael intermediate. Computation of a 2–D PMF failed to find a Diels-Alder transition structure, which was estimated to be at least 17 and 12 kcal/mol less stable than the Michael and aldol transition states. Thus, the QM/MM calculations found that the Michael-aldol mechanism is much preferred and that MPS is not a natural Diels-Alderase. For further resolution of the controversy,45–47 Hilvert and coworkers recently established that MPS efficiently promotes aldol-type chemistry lending strong support for the stepwise mechanism.48
FAAH is an integral membrane protein involved in endocannabinoid metabolism.49 The enzyme is the only characterized mammalian member of a class of serine hydrolases that bear a unique catalytic triad, Ser-Ser-Lys.50 Remarkably, FAAH hydrolyzes amides and esters with similar rates; however, the normal preference for esters reemerges when Lys142 is mutated to alanine.51 To clarify the mechanisms and unusual selectivity, the PDDG/PM2-based QM/MM approach was applied to obtain free-energy barriers for the pathways in Figure 4.52 For wild-type FAAH and oleamide, the preferred mechanism (A) was computed to involve proton transfer from Ser217 to Lys142, followed by rate-determining formation of the C–O bond between Ser241 and the substrate with concomitant proton transfer from Ser241 to Ser217. An alternative mechanism (B) that does not involve Lys142 was found to have a 4.5-kcal/mol higher barrier. Simulations were then carried out for the hydrolysis of methyl oleate; mechanism A was again preferred over B by 4.1 kcal/mol. To assess the possibility that mechanism B occurs with the Lys142Ala mutant, which lacks the basic side chain on residue 142 needed for mechanism A, hydrolysis of both substrates by the mutant was also modeled. The experimentally observed rate effects of the mutation were well-reproduced by following mechanism B. This process is a striking, multi-step sequence with simultaneous proton transfer from Ser241 to Ser217, attack of Ser241 on the carbonyl carbon of the substrate, and elimination of the leaving group and its protonation by Ser217.
In addition to the mechanistic insights, a significant technical advance was reported for the treatment of proton transfer reactions. In view of the number of such reactions for the FAAH investigation, use of traditional PMF methods would have been prohibitively resource-consuming. For a typical proton transfer, O-H…O’ → O…H–O’, it is found that the O…O’ distance remains relatively constant and that r(O-H) – r(H–O’) can be used to compute a 1-D PMF.52 This normally requires approximately 30 double-wide FEP windows using 0.02 Å Δr increments, while ca. 900 windows would be needed for a 2–D PMF using two distances as reaction coordinates. Furthermore, for proton transfers, it was found that the free-energy changes for individual windows can be fit almost perfectly by a cubic polynomial. Analytical integration yields a quartic polynomial for the overall proton-transfer PMF. It was shown that this permits accurate construction of the full PMF using only 7 FEP windows instead of the usual 30. This also allowed much more facile construction of 2-D PMFs for cases where one coordinate is a proton transfer and the other is a bond formation, e.g., for 1b → 2b in Figure 4. Nevertheless, the investigation of the reactions of FAAH entailed ca. 500 million single-point QM calculations in the course of the on-the-fly QM/MM/MC calculations.52 The need for fast, accurate QM methods is apparent.
An efficient enzyme design protocol featuring joint theoretical and experiment methodology was recently employed to create artificial enzymes for the Kemp elimination of 5-nitrobenzisoxazole.53,54 Designed enzymes were synthesized and found to yield kcat/kuncat values in the 102 – 105 range. Concurrent PDDG/PM3-based QM/MM/FEP calculations provided a deeper understanding of the enzyme structure-function relationships and gave guidance for further optimization of the catalytic performance.55 Several designed enzymes were modeled beginning from structures generated by Rosetta;54 a crystal structure for one, KE07, was subsequently obtained and found to be very close to the computational predictions. The catalytic mechanism for designs KE07, KE10, and KE15 was computed to be concerted with proton transfer generally more advanced in the transition state than breaking of the isoxazolyl N–O bond (Figure 5). The computed activation barriers near 10 kcal/mol for Kemp elimination by these enzymes were significantly smaller than the barrier of 19.8 kcal/mol computed for the reference reaction catalyzed by hydroxide ion.55 Thus, all three enzymes were anticipated to be active, as was subsequently established.
The same methodology was also applied to examine the catalysis of a Kemp elimination by a catalytic antibody, 34E4, and its Glu50Asp variant.56 The observed 30-fold rate reduction owing to this small change in the positioning of the catalytic base was well reproduced.57
Significant progress has been made in accurate QM/MM modeling of both solution-phase organic and enzymatic reactions through methodology featuring the PDDG/PM3 semiempirical QM method coupled with MC/FEP calculations. Reaction paths and activation barriers have been characterized for numerous organic reactions in multiple solvents and for several enzyme-catalyzed processes. The accord with available experimental data has been good and shows significant improvement over results from QM approaches with continuum solvent models. Nevertheless, challenges remain, especially for further improvement of fast QM methods, more rapid mapping of free-energy surfaces, and for enhanced configurational sampling of biomolecules.
Gratitude is expressed to the National Science Foundation (CHE-0446920), National Institutes of Health (GM32136), and DARPA for support of this research, to the co-workers at Auburn and Yale, and to external collaborators, especially Profs. David Baker and Kendall N. Houk.
Orlando Acevedo is a graduate of FIU and Duquesne, and was a postdoctoral associate at Yale; he is currently an Assistant Professor of Chemistry at Auburn University. Bill Jorgensen is a Sterling Professor and Director of the Division of Physical Sciences and Engineering at Yale.